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Introduction 


This  three-year  grant  was  initiated  in  December  2000.  All  work  done  under  the  auspices  of  the 
grant  is  summarized  herein,  and  the  various  papers  written  under  grant  sponsorship  are  appended 
to  provide  all  the  technical  details.  The  objective  of  this  research  was  to  provide  high-fidelity, 
finite-element-based  simulation  capability  to  assist  in  the  analysis  and  design  of  structures  for 
space  applications.  The  importance  of  such  a  research  stems  from  the  Air  Force  interest  of  space 
structures  comprising  antennas,  solar  concentrators,  sun  shields,  struts  and  booms,  and  membrane 
mirrors.  Moreover,  the  various  components  of  these  systems  often  are  in  relative  motion  with 
respect  to  one  another,  resulting  in  variable  geometry  systems  that  are  truly  multi-body  in  nature. 
Low-cost,  accurate  simulation  tools  are  essential  for  the  design,  analysis  and  assessment  of  such 
systems.  However,  efficiency  of  such  a  simulation  tool  can  only  be  achieved  through  dimensional 
reduction:  it  is  impractical  to  treat  each  structural  element  as  a  three-dimensional  (3-D)  body. 
Hence,  low-order,  accurate  models  are  needed  to  not  only  provide  the  elastic  constants  needed 
in  the  finite  element  representation  to  calculate  the  global  behavior,  but  also  to  provide  a  unified 
means  to  obtain  the  strain  and  stress  distribution  in  the  structures. 

There  are  two  main  research  aspects  in  this  project:  a  local  through-the-thickness  model  fo¬ 
cuses  on  the  evaluation  of  elastic  constants  and  displacement/strain/stress  recovery,  and  a  finite- 
element-based,  nonlinear  flexible  multi-body  dynamics  simulation  tools  based  on  those  models. 
The  basic  research  issues  that  have  been  addressed  under  the  first  aspect  are 

1.  Constructing  an  accurate  Reissner-Mindlin  type  model  for  composite  plate  using  the  varia¬ 
tional  asymptotic  method  (VAM)  [1]; 

2.  Constructing  of  an  accurate  Reissner-Mindlin  type  model  for  composite  shells  using  VAM 
including  both  geometrical  correction  due  to  initial  curvature  and  transverse  shear  effects; 

3.  Developing  a  geometrically  exact  nonlinear  shear  deformation  shell  theory  to  be  compatible 
with  the  models  developed; 

4.  Modeling  regular  composite  plates  and  shells  including  hygrothermal  effects  so  that  the 
change  of  structural  behavior  due  to  moisture  and  temperature  can  be  analyzed; 

5.  Modeling  smart  composite  plates  and  shells  made  with  piezoelectric  material. 

Approach 

The  variational  as)anptotic  method  (VAM)  has  applied  to  develop  a  high-fidelityAow-order  model 
for  composite  plates  and  shells.  Instead  of  employing  ad-hoc  assumptions,  the  VAM  relies  on  the 
use  of  the  small  parameters  that  are  inherent  to  the  structure.  This  ensures  construction  of  a  theory 
with  the  minimum  complexity  for  a  given  level  of  fidelity.  The  original  3-D  nonlinear  problem  is 
formulated  based  on  a  set  of  two-dimensional  (2-D)  intrinsic  variables  and  warping  functions  rep¬ 
resenting  the  arbitrary  deformation  of  the  normal  line.  Then  the  VAM  is  used  to  rigorously  split  the 
3-D  problem  into  two  problems:  a  nonlinear,  2-D,  plate  or  shell  analysis  over  the  reference  surface 
to  obtain  the  global  deformation  and  a  one-dimensional  (1-D)  linear  analysis  through  the  thick¬ 
ness  to  provide  the  2-D  generalized  constitutive  law  and  the  recovering  relations  to  approximate 
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the  original  3-D  results.  The  non-uniqueness  of  asymptotic  theory  correct  up  to  a  certain  order  is 
used  to  cast  the  obtained  asymptotically  correct  second-order  energy  into  a  Reissner-Mindlin  type 
model  to  accoimt  for  transverse  shear  deformation.  All  the  developed  theories  are  implemented 
in  a  finite  element  code,  Variational  Asymptotic  Plate  and  Shell  Analysis  (VAPAS).  Results  from 
VAPAS  have  been  compared  with  the  exact  solutions  available  in  the  literature,  classical  lamina¬ 
tion  theory  and  first-order  shear-deformation  theory  to  demonstrate  the  accuracy  and  power  of  the 
developed  theory. 


Work  Accomplished 

Asymptotic  Construction  of  Composite  Plate  Model 

The  development  starts  with  formulation  of  3-D  anisotropic  elasticity  problem  in  which  the  defor¬ 
mation  of  the  reference  surface  is  expressed  in  terms  of  intrinsic  2-D  variables.  The  VAM  is  then 
used  to  rigorously  split  this  3-D  problem  into  a  linear  1-D  analysis  and  a  nonlinear  2-D  “plate” 
analysis  accounting  for  transverse  shear  deformation.  The  through-the-thickness  analysis  provides 
a  constitutive  law  between  the  generalized,  two-dunensional  strains  and  stress  resultants  as  well 
as  recovering  relations  to  approximately  but  accurately  express  the  3-D  displacement,  strain  and 
stress  fields  in  terms  of  plate  variables  calculated  in  the  “plate”  analysis.  It  is  known  that  more 
than  one  theory  may  exist  that  is  asymptotically  correct  to  a  given  order.  This  non-umqueness  is 
used  to  cast  a  strain  energy  functional  &at  is  asymptotically  correct  through  the  second  order  into 
a  simple  “Reissner-Mindlin”  type  plate  theory.  Although  it  is  not  possible  in  general  to  constract 
an  asymptotically  correct  Reissner-Mindlin  type  composite  plate  theory,  an  optimization  proce¬ 
dure  is  used  to  drive  the  present  theory  as  close  to  being  asymptotically  correct  as  possible  while 
maintaining  the  simplicity  of  the  Reissner-Mindlin  formulation.  This  theory  is  firstly  formulated 
analytically  and  reported  on  an  ASME  conference  [2]  and  later  archived  in  [3].  Later,  for  the  pur¬ 
pose  to  connecting  with  2-D  finite-element  analysis  code  and  efficiency  for  multilayer  analysis,  a 
finite-element  formulation  is  developed  and  published  in  [4]. 

Asymptotic  Construction  of  Composite  Shell  Model  ^ 

A  rigorous  and  systematic  dimensional  reduction  of  a  shell-like  stmcture  is  undertaken.  Instead  of 
the  transverse  shear  refinement  we  have  for  plates,  an  accurate  shell  model  will  also  include  the 
geometry  correction  due  to  the  initial  curvatures  of  the  shell  structure.  There  are  two  main  small 
parameters,  hfl  and  h/R,  where  h  is  die  thickness  of  die  shell,  I  is  the  wavelength  of  in-plane  de¬ 
formation  and  R  is  the  minimum  radius  of  curvature  for  the  shell  structure.  Tlie  asymptotic  orders 
of  these  two  small  parameters  are  carefully  chosen  so  that  an  energy  asymptotically  correct  up  to 
the  first  order  of  h/R  and  the  second  order  of  h/l  has  been  successfully  constructed,  Our  studies 
[5, 6]  show  that  the  constructed  model  can  predict  the  stress  distribution  through  the  thickness  very 
accurately  for  shall  shells.  Improvement  to  this  theory  can  be  made  to  include  the  second  order 
of  h/J?  so  that  deep  shells  can  also  be  accurately  modeled.  A  2-D  geometrically  exact  nonlinear 
shear  deformation  shell  theory  [7]  has  been  developed  to  be  consistent  with  the  Reissner-Mindlin 
type  model  constructed.  A  complete  set  of  kinematical  and  intrinsic  equilibrium  equations  are 
derived  for  shells  undergoing  large  displacements  and  rotations  but  with  small,  2-D  generalized 


3 


strains.  The  large  rotation  is  represented  by  the  general  finite  rotation  is  along  the  normal  line.  It  is 
shown  that  the  rotation  of  the  frame  about  the  normal  line  is  not  zero  and  that  it  can  be  expressed  in 
terms  of  other  global  deformation  variables.  It  is  also  shown  that  only  five  equilibriiun  equations 
can  be  derived  in  this  manner  because  the  component  of  virtual  rotation  about  the  normal  is  not 
independent. 

A  Thermolelastic  Composite  Plate  Model 

Composite  structures  are  more  sensitive  and  vulnerable  to  temperature  change  than  their  isotropic 
coimterpart  because  the  thermal  expansion  coefficients  of  different  constituents  of  the  material  are 
usually  dramatically  different  from  each  other  resulting  in  high  stresses  due  to  sudden  temperature 
change.  The  analysis  including  thermal  effects  is  much  more  involved  than  that  for  isothermal 
conditions.  A  Reissner-Mindlin  type  plate  model  capable  of  performing  a  thermoelastic  stress 
analysis  of  laminated  composite  plates  has  been  constructed  by  the  VAM.  It.  is  shown  in  [8]  that 
although  the  resulting  theory  is  of  the  simple  Reissner-Mindlin  form,  it  has  an  accuracy  comparable 
to  a  higher-order  layerwise  theory.  The  hygro  effect  due  to  moisture  to  composite  plates  can  also 
handled  in  exactly  the  same  procedure  except  one  has  to  replace  the  thermal  expansion  coefficients 
with  hygroscopic  expansion  coefficients  and  temperature  with  moisture. 

A  Thermopiezoelastic  Composite  Plate  Model 

A  Reissner-Mindlin  type  model  for  analyzing  laminated  smart  composite  plates  including  piezo¬ 
electric  layers  under  mechanical,  thermal  and  electric  loads  has  been  developed.  This  model  can 
analyze  the  one-way  coupling  between  stracture  and  thermal,  electric  field.  The  non-mechanical 
stress  resultants  due  to  temperature  and  actuation  can  be  predicted  using  this  model.  This  theory 
has  been  implemented  into  the  computer  program  VAPAS  and  validated  against  published  exact 
solutions  for  simple  problems.  As  reported  in  [9, 10],  the  recovered  3-D  stress  distribution  due  to 
temperature,  or  electricity  has  a  good  agreement  with  the  exact  solutions. 

A  Thermopiezoelastic  Composite  Shell  Model 

Finally,  the  above  developed  thermopiezoelastic  model  for  smart  plates  is  extended  to  model  smart 
shells  made  embedded  with  piezoelectric  layers.  All  the  effects  due  to  temperature,  electricity, 
initial  curvature  and  transverse  shear  are  represented  in  this  model.  This  model  is  also  implemented 
in  the  computer  program  VAPAS  and  a  paper  giving  the  results  has  been  accepted  for  presentation 
at  the  2004  AIAA  SDM  conference  [1 1]. 

Multibody  Modeling  of  Shells  and  Membranes 

The  classical  approach  to  the  numerical  simulation  of  flexible  multibody  systems  proceeds  in  two 
steps;  first,  the  equations  of  motion  of  the  system  are  written  in  a  convenient  form,  then  general 
purpose  Differential  Algebraic  Equations  (DAF.)  solvers  are  used  tn  integrate  rhest*  Rgn^tinT^s  in 
the  time  domain.  General  purpose  DAE  integrators  are  specifically  designed  for  effectively  dealing 
with  the  dual  differential/algebraic  nature  of  the  equations  but  are  otherwise  unaware  of  the  specific 
features  and  characteristics  of  the  equations  being  solved. 
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The  equations  governing  nonlinear  flexible  multibody  systems  with  shells  and  membranes 
present  very  specific  features.  First,  they  are  characterized  by  linear  and  rotational  tensorial  fields 
describing  kinematic  (displacements,  velocities)  and  co-kinematic  (forces,  momenta)  quantities. 
Second,  nonlinearities  can  arise  from  several  sources:  large  displacements  and  finite  rotations 
(geometric  nonlinearities),  or  nonlinear  constitutive  laws  for  the  deformable  components  of  the 
system  (material  nonlinearities).  Third,  a  distinguishing  feature  of  multibody  systems  is  the  pres¬ 
ence  of  joints  which  impose  different  types  of  kinematic  constraints  between  the  various  bodies  of 
the  system.  More  often  than  not,  constraints  are  modeled  via  the  Lagrange  multipliers  technique 
that  imposes  the  nonlinear  algebraic  constraints  on  the  system.  Fourth,  the  exact  solution  of  the 
equations  of  motion  implies  the  preservation  of  a  number  of  dynamic  invariants,  such  as  energy 
and  momenta.  Fifth,  when  the  elastic  bodies  of  the  system  are  modeled  by  means  of  an  appro¬ 
priate  spatial  discretization  process,  such  as  the  finite  element  method,  high  fi'equency  modes  are 
introduced  in  the  system.  Finally,  when  dealing  with  shells  and  membranes,  the  ratio  of  in-plane 
to  out-of-plane  fi'equencies  is  extremely  high,  leading  to  very  stiff  systems. 

While  standard  approaches  perform  adequately  for  a  number  of  simulations,  problems  can  arise 
when  modeling  flexible,  nonlinear  multibody  systems  involving  shells  and  membranes.  In  this 
case,  robust  algorithms  that  satisfy  precise  requirements  should  be  designed  fox  the  time  integration 
of  such  systems.  A  formulation  of  time  inte^tors  for  shells  and  membranes  was  developed  that 
satisfies  the  following  requirements:  nonlinear  unconditional  stability  of  the  scheme,  a  rigorous 
treatment  of  all  nonlinearities,  the  exact  satisfaction  of  the  constraints,  and  the  presence  of  high 
frequency  numerical  dissipation.  The  proof  of  nonlinear  unconditional  stability  stems  from  two 
physical  characteristics  of  multibody  systems  that  will  be  reflected  in  the  numerical  scheme:  the 
preservation  of  the  total  mechanical  energy,  and  the  vanishing  of  the  work  performed  by  constraint 
forces.  Numerical  dissipation  is  obtained  by  letting  the  solution  drift  from  the  constant  energy 
manifold  in  a  controlled  manner  in  such  a  way  that  at  each  time  step,  energy  can  be  dissipated  but 
not  created. 

A  novel  integration  scheme  for  nonlinear  dynamics  of  geometrically  exact  shells  was  devel¬ 
oped  based  on  the  inextensible  director  assumption.  The  new  algorithm  is  designed  so  as  to  imply 
the  strict  decay  of  the  system  total  mechanical  energy  at  each  time  step,  and  consequently  imcondi- 
tional  stability  is  achieved  in  the  nonlinear  regime.  Furthermore,  the  scheme  features  tunable  high 
frequency  munerical  damping  and  it  is  therefore  stiflfiy  accurate.  The  method  is  tested  for  a  finite 
element  spatial  formulation  of  shells  based  on  mixed  interpolations  of  strain  tensorial  components 
and  on  a  two-parameter  representation  of  director  rotations.  The  robustness  of  the  scheme  is  illus¬ 
trated  with  the  help  of  numerical  examples  presented  at  conferences  [12,  13]  and  later  published 
in  a  paper  [14]. 

Next,  these  algorithms  were  extended  to  deal  with  the  complexities  associated  with  multibody 
systems.  Results  were  first  presented  at  conferences  [15,16]  and  later  published  in  papers  [17,18] 
addressing  various  aspects  of  the  problem.  The  modeling  of  cables  and  membranes  was  addressed 
in  [19]. 


Conclusions  and  Suggestions 

A  rigorous  framework  to  constract  accurate  model  for  composite  plates  and  shells  has  been  estab¬ 
lished  under  this  grant  using  the  variational  asymptotic  method.  Models  constracted  this  way  can 
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not  only  provide  accurate  2-D  constitutive  model  for  plate  or  shell  analysis  but  also  a  consistent 
way  to  recover  the  original  3-D  displacement/strain/stress  fields.  All  the  theories  developed  have 
been  implemented  in  the  computer  program  VAPAS.  VAPAS  has  been  connected  with  DYMORE 
to  provide  an  efficient  yet  accurate  simulation  for  complex  aerospace  systems. 

More  work  needs  to  be  done  to  verify  the  developed  and  implemented  thermopiezoelastic 
model  for  smart  shells  against  available  published  results.  And  also  the  developed  shell  mod¬ 
els  could  be  improved  by  constructing  an  energy  asymptotically  correct  up  to  the  second  order  in 
h/K  However,  this  requires  significant  research  effort  and  would  also  require  the  development  of 
shell  models  capable  of  providing  the  required  inputs  to  such  models. 
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Summary 

High-fidelity,  finite-element-based  dimensionally  reduced  models  have  been  constructed  for  com¬ 
posite  plates  and  shells,  including  hygrothermal  and  piezoelectric  effects,  under  the  sponsorship 
of  AFOSR.  In  these  models,  the  smallness  of  the  thickness  has  been  used  to  advantage  to  rig¬ 
orously  reduce  the  original  three-dimensional  geometrically  nonlinear  elasticity  theory  to  two- 
dimensional  Reissner-Mindlin  type  theory  for  plates  and  shells.  The  resulting  theory  can  achieve 
an  accuracy  comparable  to  higher-order  layerwise  theories  at  the  cost  of  only  a  first-order  shear 
deformation  theory.  The  dimensional  reduction  process  and  the  recovery  relations  for  the  original 
three-dimensional  displacements/strains/stresses  are  implemented  in  a  finite-element  code,  Varia¬ 
tional  Asymptotic  Plate  and  Shell  Analysis  (VAPAS).  This  program  is  connected  with  DYMORE, 
a  nonlinear  finite-element  based  multi-body  dynamic  code  to  provide  an  efficient  and  accurate  sim¬ 
ulation  capability  for  space  systems  involving  composite  and  inflatable  components  with  actuated 
elements  which  are  required  for  current  and  future  Air  Force  missions. 


Introduction 


This  three-year  grant  was  initiated  in  December  2000.  All  work  done  under  the  auspices  of  the 
grant  is  summarized  herein,  and  the  various  papers  written  under  grant  sponsorship  are  appended 
to  provide  all  the  technical  details.  The  objective  of  this  research  was  to  provide  high-fidelity, 
finite-element-based  simulation  capability  to  assist  in  the  analysis  and  design  of  structures  for 
space  applications.  The  importance  of  such  a  research  stems  from  the  Air  Force  interest  of  space 
structures  comprising  antennas,  solar  concentrators,  sun  shields,  struts  and  booms,  and  membrane 
mirrors.  Moreover,  the  various  components  of  these  systems  often  are  in  relative  niotion  with 
respect  to  one  another,  resulting  in  variable  geometry  systems  that  are  truly  multi-body  in  nature. 
Low-cost,  accurate  simulation  tools  are  essential  for  the  design,  analysis  and  assessment  of  such 
systems.  However,  efficiency  of  such  a  simulation  tool  can  only  be  achieved  through  dimensional 
reduction;  it  is  impractical  to  treat  each  structural  element  as  a  three-dimensional  (3-D)  body. 
Hence,  low-order,  accurate  models  are  needed  to  not  only  provide  the  elastic  constants  needed 
in  the  finite  element  representation  to  calculate  the  global  behavior,  but  also  to  provide  a  unified 
means  to  obtain  the  strain  and  stress  distribution  in  the  structures. 

There  are  two  main  research  aspects  in  this  project:  a  local  through-the-thickness  model  fo¬ 
cuses  on  the  evaluation  of  elastic  constants  and  displacement/strain/stress  recovery,  and  a  finite- 
element-based,  nonlinear  flexible  multi-body  dynamics  simulation  tools  based  on  those  models. 
The  basic  research  issues  that  have  been  addressed  under  the  first  aspect  are 

1.  Constructing  an  accurate  Reissner-Mindlin  type  model  for  composite  plate  using  the  varia¬ 
tional  asymptotic  method  (VAM)  [IJ; 

2.  Constructing  of  an  accurate  Reissner-Mindlin  type  model  for  composite  shells  using  VAM 
including  both  geometrical  correction  due  to  initial  curvature  and  transverse  shear  effects; 

3.  Developing  a  geometrically  exact  nonlinear  shear  deformation  shell  theory  to  be  compatible 
with  the  models  developed; 

4.  Modeling  regular  composite  plates  and  shells  including  hygrothermal  effects  so  that  the 
change  of  structural  behavior  due  to  moisture  and  temperature  can  be  analyzed; 

5.  Modeling  smart  composite  plates  and  shells  made  with  piezoelectric  material. 

Approach 

The  variational  asymptotic  method  (VAM)  has  applied  to  develop  a  high-fidelity/low-order  model 
for  composite  plates  and  shells.  Instead  of  employing  ad-hoc  assumptions,  the  VAM  relies  on  the 
use  of  the  small  parameters  that  are  inherent  to  the  structure.  This  ensures  construction  of  a  theory 
with  the  minimum  complexity  for  a  given  level  of  fidelity.  The  original  3-D  nonlinear  problem  is 
formulated  based  on  a  set  of  two-dimensional  (2-D)  intrinsic  variables  and  warping  functions  rep¬ 
resenting:  the  arbitrary  deformation  of  the  normal  line.  Then  the  VAM  is  used  to  rigorously  split  the 
3-D  problem  into  two  problems:  a  nonlinear,  2-D,  plate  or  shell  analysis  over  the  reference  surface 
to  obtain  the  global  deformation  and  a  one-dimensional  (i-D)  linear  analysis  through  the  thick¬ 
ness  to  provide  the  2-D  generalized  constitutive  law  and  the  recovering  relations  to  approximate 


the  original  3-D  results.  The  non-uniqueness  of  asymptotic  theory  correct  up  to  a  certain  order  is 
used  to  cast  the  obtained  asymptotically  correct  second-order  energy  into  a  Reissner-Mindlin  type 
model  to  account  for  transverse  shear  deformation.  All  the  developed  theories  are  implemented 
in  a  finite  element  code.  Variational  Asymptotic  Plate  and  Shell  Analysis  (VAPAS).  Results  from 
VAPAS  have  been  compared  with  the  exact  solutions  available  in  the  literature,  classical  lamina¬ 
tion  theory  and  first-order  shear-deformation  theory  to  demonstrate  the  accuracy  and  power  of  the 
developed  theory. 

Work  Accomplished 

Asymptotic  Construction  of  Composite  Plate  Model 

The  development  starts  with  formulation  of  3-D  anisotropic  elasticity  problem  in  which  the  defor¬ 
mation  of  the  reference  surface  is  expressed  in  terms  of  intrinsic  2-D  variables.  The  YAM  is  then 
used  to  rigorously  split  this  3-D  problem  into  a  linear  1-D  analysis  and  a  nonlinear  2-D  “plate” 
analysis  accounting  for  transverse  shear  deformation.  The  through-the-thickness  analysis  provides 
a  constitutive  law  between  the  generalized,  two-dimensional  strains  and  stress  resultants  as  well 
as  recovering  relations  to  approximately  but  accurately  express  the  3-D  displacement,  strain  and 
stress  fields  in  terms  of  plate  variables  calculated  in  the  “plate”  analysis.  It  is  known  that  more 
than  one  theory  may  exist  that  is  asymptotically  correct  to  a  given  order.  This  non-uniqueness  is 
used  to  cast  a  strain  energy  functional  that  is  asymptotically  correct  through  the  second  order  into 
a  simple  “Reissner-Mindlin”  type  plate  theory.  Although  it  is  not  possible  in  general  to  construct 
an  asymptotically  correct  Reissner-Mindlin  type  composite  plate  theory,  an  optimization  proce¬ 
dure  is  used  to  drive  the  present  theory  as  close  to  being  asyrriptotically  correct  as  possible  while 
maintaining  the  simplicity  of  the  Reissner-Mindlin  formulation.  This  theory  is  firstly  formulated 
analytically  and  reported  on  an  ASME  conference  [2J  and  later  archived  in  13|.  Later,  for  the  pur¬ 
pose  to  connecting  with  2-D  finite-element  analysis  code  and  efficiency  for  multilayer*analysis,  a 
finite-element  formulation  is  developed  and  published  in  [4]. 

Asymptotic  Construction  of  Composite  Shell  Model 

A  rigorous  and  systematic  dimensional  reduction  of  a  shell-like  structure  is  undertaken.  Instead  of 
the  transverse  shear  refinement  we  have  for  plates,  an  accurate  shell  model  will  also  include  the 
geometry  correction  due  to  the  initial  curvatures  Of  the  shell  structure.  There  are  two  main  small 
parameters,  h/l  and  h/R,  where  h  is  the  thickness  of  the  shell,  I  is  the  wavelength  of  in-plane  de¬ 
formation  and  Ris  the  minimum  radius  of  curvature  for  the  shell  structure.  The  asymptotic  orders 
of  these  two  small  parameters  are  carefully  chosen  so  that  an  energy  asymptotically  correct  up  to 
the  first  order  of  h/R  and  the  second  order  of  h/l  has  been  successfully  constructed.  Our  studies 
(5, 6|  show  that  the  constructed  model  can  predict' the  stress  distribution  through  the  thickness  very 
accurately  for  shall  shells.  Improvement  to  this  theory  can  be  made  to  include  the  second  order 
of  h/R  so  that  deep  shells  can  also  be  accurately  modeled.  A  2-D  geometrically  exact  nonlinear 
shear  deformation  shell  theory  [71  has  been  developed  to  be  consistent  with  the  Reissner-Mindlin 
type  model  constructed.  A  complete  set  of  kinematical  and  intrinsic  equilibrium  equations  are 
derived  for  shells  undergoing  large  displacements  and  rotations  but  with  small,  2-D  generalized 


strains.  The  large  rotation  is  represented  by  the  general  finite  rotation  is  along  the  normal  line.  It  is 
shown  that  the  rotation  of  the  frame  about  the  normal  line  is  not  zero  and  that  it  can  be  expressed  in 
terms  of  other  global  deformation  variables.  It  is  also  shown  that  only  five  equilibrium  equations 
can  be  derived  in  this  manner  because  the  component  of  virtual  rotation  about  the  normal  is  not 
independent. 

A  Thermolelastic  Composite  Plate  Model 

Composite  structures  are  more  sensitive  and  vulnerable  to  temperature  change  than  their  isotropic 
counterpart  because  the  thermal  expansion  coefficients  of  different  constituents  of  the  material  are 
usually  dramatically  different  from  each  other  resulting  in  high  stresses  due  to  sudden  temperature 
change.  The  analysis  including  thermal  effects  is  much  more  involved  than  that  for  isothermal 
conditions.  A  Reissner-Mindlin  type  plate  model  capable  of  performing  a  thermoelastic  stress 
analysis  of  laminated  composite  plates  has  been  constructed  by  the  VAM.  It  is  shown  in  [8]  that 
although  the  resulting  theory  is  of  the  simple  Reissner-Mindlin  form,  it  has  an  accuracy  comparable 
to  a  higher-order  layerwise  theory.  The  hygro  effect  due  to  moisture  to  composite  plates  can  also 
handled  in  exactly  the  same  procedure  except  one  has  to  replace  the  thermal  expansion  coefficients 
with  hygroscopic  expansion  coefficients  and  temperature  with  moisture. 


A  Thermopiezoellastic  Composite  Plate  Model 

A  Reissner-Mindlin  type  model  for  analyzing  laminated  smart  composite  plates  including  piezo¬ 
electric  layers  under  mechanical,  thermal  and  electric  loads  has  been  developed.  This  model  can 
analyze  the  one-way  coupling  between  structure  and  thermal,  electric  field.  The  non-mechanical 
stress  resultants  due  to  temperature  and  actuation  can  be  predicted  using  this  model.  This  theory 
has  been  implemented  into  the  computer  program  VAPAS  and  validated  against  published  exact 
solutions  for  simple  problems.  As  reported  in  [9,  iOJ,  the  recovered  3-D  stress  distribution  due  to 
temperature,  br  electricity  has  a  good  agreement  with  the  exact  solutions. 

A  Thermopiezoelastic  Composite  Shell  Model 

Finally,  the  above  developed  thermopiezoelastic  model  for  smart  plates  is  extended  to  model  smart 
shells  made  embedded  with  piezoelectric  layers.  All  the  effects  due  to  temperature,  electricity, 
initial  curvature  and  transverse  shear  are  represented  in  this  model.  This  model  is  also  implemented 
in  the  computer  program  VAPAS  and  a  paper  giving  the  results  has  been  accepted  for  presentation 
at  the  2004  AIAA  SDM  conference  [1 1 1. 

Multibody  Modeling  of  Shells  and  Membranes 

The  classical  approach  to  the  numerical  simulation  of  flexible  multibody  systems  proceeds  in  two 
steps:  first,  the  equations  of  motion  of  the  system  are  written  in  a  convenient  form,  then  general 
purpose  Differential  Algebraic  Equations  (DAE)  solvers  are  used  to  integrate  these  equations  in 
the  time  domain.  General  purpose  DAE  integrators  are  specifically  designed  for  effectively  dealing 
with  the  dual  differential/algebraic  nature  of  the  equations  but  are  otherwise  unaware  of  the  specific 
features  and  characteristics  of  the  equations  being  solved. 


The  equations  governing  nonlinear  flexible  multibody  systems  with  shells  and  membranes 
present  very  specific  features.  First,  they  are  characterized  by  linear  and  rotational  tensorial  fields 
describing  kinematic  (displacements,  velocities)  and  co-kinematic  (forces,  momenta)  quantities. 
Second,  nonlinearities  can  arise  from  several  sources:  large  displacements  and  finite  rotations 
(geometric  nonlinearities),  or  nonlinear  constitutive  laws  for  the  deformable  components  of  the 
system  (material  nonlinearities).  Third,  a  distinguishing  feature  of  multibody  systems  is  the  pres¬ 
ence  of  joints  which  impose  different  types  of  kinematic  constraints  between  the  various  bodies  of 
the  system.  More  often  than  not,  constraints  are  modeled  via  the  Lagrange  multipliers  technique 
that  imposes  the  nonlinear  algebraic  constraints  on  the  system.  Fourth,  the  exact  solution  of  the 
equations  of  motion  implies  the  preservation  of  a  nunriber  of  dynamic  invariants,  such  as  energy 
and  momenta.  Fifth,  when  the  elastic  bodies  of  the  system  are  modeled  by  means  of  an  appro¬ 
priate  spatial  discretization  process,  such  as  the  finite  element  method,  high  frequency  modes  are 
introduced  in  the  system.  Finally,  when  dealing  with  shells  and  membranes,  the  ratio  of  in-plane 
to  out-of-plane  frequencies  is  extremely  high,  leading  to  very  stiff  systems. 

While  standard  approaches  perform  adequately  for  a  number  of  simulations,  problems  can  arise 
when  modeling  flexible,  nonlinear  multibody  systems  involving  shells  and  membranes.  In  this 
case,  robust  algorithms  that  satisfy  precise  requirements  should  be  designed  for  the  time  integration 
of  such  systems.  A  formulation  of  time  integrators  for  shells  and  membranes  was  developed  that 
satisfies  the  following  requirements:  nonlinear  unconditional  stability  of  the  scheme,  a  rigorous 
treatment  of  all  nonlinearities,  the  exact  satisfaction  of  the  constraints,  and  the  presence  of  high 
frequency  numerical  dissipation.  The  proof  of  nonlinear  unconditional  stability  stems  from  two 
physical  characteristics  of  multibody  systems  that  will  be  reflected  in  the  numerical  scheme:  the 
preservation  of  the  total  mechanical  energy,  and  the  vanishing  of  the  work  performed  by  constraint 
forces.  Numerical  dissipation  is  obtained  by  letting  the  solution  drift  from  the  constant  energy 
manifold  in  a  controlled  manner  in  such  a  way  that  at  each  time  step,  energy  can  be  dissipated  but 
not  created. 

A  novel  integration  scheme  for  nonlinear  dynamics  of  geometrically  exact  shells  was  devel¬ 
oped  based  on  the  inextensible  director  assumption.  The  new  algorithm  is  designed  so  as  to  imply 
the  strict  decay  of  the  system  total  mechanical  energy  at  each  time  step,  and  consequently  uncondi¬ 
tional  stability  is  achieved  in  the  nonlinear  regime.  Furthermore,  the  scheme  features  tunable  high 
frequency  numerical  damping  and  it  is  therefore  stiffly  accurate.  The  method  is  tested  for  a  finite 
elerhent  spatial  formulation  of  shells  based  on  mixed  interpolations  of  strain  tensorial  components 
and  on  a  two-parameter  representation  of  director  rotations.  The  robustness  of  the  scheme  is  illus¬ 
trated  with  the  help  of  numerical  examples  presented  at  conferences  [12,  13 1  and  later  published 
in  a  paper  [14]. 

Next,  these  algorithms  were  extended  to  deal  with  the  complexities  associated  with  multibody 
systems.  Results  were  first  presented  at  conferences  [15,  16]  and  later  published  in  papers  [17,  18] 
addressing  various  aspects  of  the  problem.  The  modeling  of  cables  and  membranes  was  addressed 
in|19|. 

Conclusions  and  Suggestions 

A  rigorous  framework  to  construct  accurate  model  for  composite  plates  and  shells  has  been  estab¬ 
lished  under  this  grant  using  the  variational  asymptotic  method.  Models  constructed  this  way  can 


not  only  provide  accurate  2-D  constitutive  model  for  plate  or  shell  analysis  but  also  a  consistent 
way  to  recover  the  original  3-D  displacement/strain/stress  fields.  All  the  theories  developed  have 
been  implemented  in  the  computer  program  VAPAS.  VAPAS  has  been  connected  with  DYMORE 
to  provide  an  efficient  yet  accurate  simulation  for  complex  aerospace  systems. 

More  work  needs  to  be  done  to  verify  the  developed  and  implemented  thermopiezoelastic 
model  for  smart  shells  against  available  published  results.  And  also  the  developed  shell  mod¬ 
els  could  be  improved  by  constructing  an  energy  asymptotically  correct  up  to  the  second  order  in 
h/R.  However,  this  requires  significant  research  effort  and  would  also  require  the  development  of 
shell  models  capable  of  providing  the  required  inputs  to  such  models. 
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Appendix  D 


A  Geometrically  Nonlinear  Shear 
Deformation  Theory  for  Composite  Shells 

Wenbin  Yu*  and  Dewey  H.  Hodges^ 

Georgia  Institute  of  Technology,  Atlanta,  Georgia  30332-0150 

A  geometrically  nonlinear  shear  deformation  theory  has  been  developed  for  elastic  shells  to  ac¬ 
commodate  a  constitutive  model  suitable  for  composite  shells  when  modeled  as  a  two-dimensional 
continuum.  A  complete  set  of  kinematical  and  intrinsic  equilibrium  equations  are  derived  for  shells 
undergoing  large  displacements  and  rotations  but  with  small,  two-dimensional,  generalized  Strains. 
The  large  rotation  is  represented  by  the  general  finite  rotation  of  a  frame  embedded  in  the  unde¬ 
formed  configuration,  of  which  one  axis  is  along  the  normal  line.  The  unit  vector  along  the  normal 
line  of  the  undeformed  reference  surface  is  not  in  general  normal  to  the  deformed  reference  sur¬ 
face  because  of  transverse  shear.  It  is  shown  that  the  rotation  of  the  frame  about  the  normal  line 
is  not  zero  and  that  it  can  be  expressed  in  terms  of  other  global  deformation  variables.  Based 
on  a  generalized  constitutive  model  obtained  from  an  asymptotic  dimensional  reduction  from  the 
three-dimensional  energy,  and  in  the  form  of  a  Reissner-Mindlin  type  theory,  a  set  of  intrinsic  equi¬ 
librium  equations  and  boundary  conditions  follow.  It  is  shown  that  only  five  equilibrium  equations 
can  be  derived  in  this  manner  because  the  corhponent  of  virtual  rotation  about  the  normal  is  not 
independent.  It  is  shown,  however,  that  these  equilibrium  equations  contain  terms  that  cannot  be 
obtained  without  the  use  of  all  three  components  of  the  finite  rotation  vector. 

Introduction 

For  an  elastic  three-dimensional  (3-D)  continuum,  there  are  two  types  of  nonlinearity-  geomet¬ 
rical  and  physical.  A  theory  is  geometrically  nonlinear  if  the  kinematical  (strain-displacement) 
relations  are  nonlinear  but  the  constitutive  (stress-strain)  relations  are  linear.  This  kind  of  theory 
allows  large  displacements  and  rotations  with  the  restriction  that  strain  must  be  small.  A  physi¬ 
cally  (or  materially)  nonlinear  theory  is  necessary  for  biological,  rubber-like  or  inflatable  structures 
where  the  strain  cannot  be  considered  small,  and  a  nonlinear  constitutive  law  is  needed  to  relate 
the  stress  and  strain.  Although  this  classification  seems  obvious  and  clear  for  a  structure  modeled 
as  a  3-D  continuum,  it  becomes  somewhat  ambiguous  to  model  dimensionally  reducible  structures 
-  structures  have  one  or  two  dimensions  much  smaller  than  the  other(s)  such  as  beams,  plates  and 
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shells  [1]  -  using  reduced  one-dimensional  (1-D)  or  two-dimensional  (2-D)  models.  A  nonlin¬ 
ear  constitutive  law  for  the  reduced  structural  model  can  in  some  circumstances  be  obtained  from 
the  reduction  of  a  geometrically  nonlinear  3-D  theory.  For  example,  in  the  so-called  Wagner  or 
trapeze  effect  [2-5],  the  effective  torsional  rigidity  is  increased  due  to  axial  force.  This  physically 
nonlinear  1-D  model  stems  from  a  purely  geometrically  nonlinear  theory  at  the  3-D  level.  On  the 
other  hand,  the  present  paper  focuses  on  a  geometrically  nonlinear  analysis  at  the  3-D  level  which 
becomes  a  geometrically  nonlinear  analysis  at  the  two-dimensional  2-D  as  well.  That  is,  the  2- 
D  generalized  strain-displacement  relations  are  nonlinear  while  the  2-D  generalized  stress-strain 
relations  turn  out  to  be  linear. 

A  shell  is  a  3^D  body  with  a  relatively  small  thickness  and  a  smooth  reference  surface.  The 
feature  of  the  small  thickness  attracts  researchers  to  simplify  their  analyses  by  reducing  the  original 
3-D  problem  to  a  2-D  problem  by  taking  advantage  of  the  thinness.  By  comparison  with  the 
original  3-D  problem,  an  exact  shell  theory  does  not  exist.  Dimensional  reduction  is  an  inherently 
approximate  process.  Shell  theory  is  a  very  old  subject,  since  the  vibration  of  a  bell  was  attempted 
by  Euler  even  before  elasticity  theory  was  well  established  [6].  Even  so,  shell  theory  still  receives 
a  lot  of  attention  from  modern  researchers  because  it  is  used  so  extensively  in  so  many  engineering 
applications.  Moreover,  many  shells  are  now  made  with  advanced  materials  that  have  only  recently 
become  available. 

Generally  speaking,  shell  theories  can  be  classified  according  to  direct,  derived  and  mixed  ap¬ 
proaches.  The  direct  approach,  which  originated  with  the  Cosserat  brothers  [7],  models  a  shell 
directly  as  a  2-D  “orientated”  continuum.  Naghdi  [8]  provided  an  extensive  review  of  this  kind  of 
approach.  Although  the  direct  approach  is  elegant  and  able  to  account  for  transverse  and  normal 
strains  and  rotations  associated  with  couple  stresses,  it  nowhere  connects  with  the  fact  that  a  shell 
is  a  3-D  body  and  thus  completely  isolates  itself  from  3-D  continuum  mechanics.  This  could  be 
the  main  reason  that  this  approach  has  not  been  much  appreciated  in  the  engineering  community. 
One  of  the  complaints  of  these  approaches  that  they  are  difficult  for  numerical  implementation 
has  been  answered  by  Simo  and  his  co-workers  by  providing  an  efficient  formulation  “free  from 
mathematical  complexities  and  suitable  for  large  scale  computation”  [9, 10].  And  mOre  recently 
a  similar  theory  was  developed  by  Ibrahimbegovic  [11]  to  include  drilling  rotations  so  that  not- 
so-smooth  shell  structures  can  be  analyzed  conveniently.  However,  the  main  complaint  remains 
that  these  approaches  lack  a  meaningful  way  to  find  the  constitutive  models  “which  Can  only  be 
experienced  and  formulated  properly  in  our  3-D  real  world”  [12].  Reissner  [13]  developed  a  very 
general  nonlinear  shell  theory  introducing  twelve  generalized  strains  by  considering  the  dynamics 
of  stress  resultants  and  couples  on  the  reference  surface  as  the  basis.  He  gracefully  avoided  the 
awkwardness  of  finding  a  proper  constitutive  model  by  pointing  out  two  possible  means  to  estab¬ 
lish  them.  It  is  recommended  in  [13]  that  one  could  either  design  experiments  to  determine  the 
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constitutive  constants  without  explicit  reference  to  the  3-D  nature  of  the  structure  or  derive  an  ap¬ 
propriate  2-D  model  from  the  given  knowledge  of  the  constitutive  relations  for  the  real  3-D  model 
of  the  structure. 

Derived  approaches  reduce  the  original  3-D  elasticity  problem  into  a  2-D  problem  to  be  solved 
over  the  reference  surface.  Such  reductions  are  usually  carried  out  in  one  of  two  ways.  The  most 
common  approach  is  to  assume  a  priori  the  distribution  of  3-D  quantities  through  the  thickness  and 
then  to  construct  a  2-D  strain  energy  per  unit  area  by  integrating  the  3-D  energy  per  unit  volume 
through  the  thickness.  Remarkably,  classical  (also  known  as  Kirchhoff-Love  type  theory),  first- 
order  shear  deformation  (also  known  as  Reissner-Mindlin  type  theory),  higher-order,  and  layer- 
wise  shell  theories  all  fall  into  this  category,  including  the  theories  proposed  by  Reddy  [14],  for 
example.  Another  approach  is  to  apply  an  asymptotic  method  to  expand  all  quantities  into  an 
asymptotic  series  of  the  thickness  coordinate,  so  that  a  sequence  of  2-D  problems  can  be  solved 
according  to  the  different  orders. 

The  mixed  approach  is  used  in  [15]  based  on  the  argument  that  all  the  3-D  elasticity  equations 
except  the  constitutive  relations  are  independent  of  the  material  properties,  such  as  the  kinematical 
relations,  equilibrium  of  momentum  and  forces.  The  constitutive  law  must  be  determined  experi¬ 
mentally,  and  hence  it  is  avoidable  that  it  is  approximate.  Libai  and  Simmonds  [15]  obtain  exact 
shell  equations  for  the  balance  of  momentum,  heat  flow  and  an  entropy  inequality  from  the  3-D 
continuum  mechanics  via  integration  through  the  thickness.  An  analogous  2-D  constitutive  law  is 
postulated  due  to  the  fact  that  even  3-D  constitutive  laws  are  inexact. 

There  is  a  sense  in  which  the  present  approach  can  also  be  considered  as  mixed.  The  2-D 
constitutive  model  is  obtained  by  the  Variational  Asymptotic  Method  (VAM)  [16]  such  that  the 
2-D  energy  is  as  close  to  an  asymptotic  approximation  of  the  original  3-D  energy  as  possible 
[17].  The  process  of  constructing  the  constitutive  model  defines  the  reference  surface  and  the 
kinematics  of  this  surface  are  geometrically  exact  formulated  in  an  intrinsic  format.  The  2-D 
equilibrium  equations  are  obtained  from  the  2-D  energy  with  the  knowledge  of  the  variations  of 
the  generalized  strains.  The  only  approximate  part  of  our  2-D  shell  theory  is  the  constitutive  law 
which  is  not  postulated  but  is  mathematically  obtained  by  VAM. 

Shell  Kinematics 

The  equations  of  2-D  shell  theory  are  written  over  the  domain  of  the  reference  surface,  on 
which  every  point  can  be  represented  by  a  position  vector  r  in  the  undeformed  configuration  and 
R  in  the  deformed  configuration  (see  Fig.  1)  with  respect  to  a  fixed  point  O  in  the  space.  A  set 
of  two  curvilinear  coordinates,  Xa,  are  required  to  located  a  point  on  the  reference  surface.  The 
coordinates  are  so-called  convected  coordinates  such  that  every  point  of  the  configuration  has  the 
same  coordinates  during  the  deformation.  (Here  and  throughout  the  paper  Latin  indices  assume  1, 
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2. 3;  and  Greek  indices  assume  values  1  and  2.  Dummy  indices  are  summed  over  their  range  except 
where  explicitly  indicated.)  Without  loss  of  generality,  Xq  are  chosen  to  be  the  lines  of  curvatures 
of  the  surface  to  simplify  the  formulation.  For  the  purpose  of  representing  finite  rotations,  an 
orthonormal  triad  bj  is  introduced  for  the  initial  configuration,  such  that 

ha  =  aia/Aa  bs  =  bi  X  b2  (1) 

where  is  the  set  of  base  vectors  associated  with  Xa  and  ..4^  are  the  Lame  parameters,  defined  as 

&a  ~  ^,a  Aa  —  \/^a  '  (2) 


From  the  differential  geometry  of  the  surface  and  following  [13]  and  [18]  one  can  express  the 
derivatives  of  b,  as 

bj_£j  =  .4£ikQ  X  bj  (3) 

where  is  the  curvature  vector  measured  in  bj  with  the  components 


fCa  —  |_  fCa2  ^asj 


(4) 


in  which  refers  to  out-of-plane  curvatures.  We  note  that  A’la  =  A21  =  0  because  the  coordinates 
are  the  tines  of  curvatures.  The  geodesic  curvatures  kas  can  be  e.xpressed  in  terms  of  the  Lame 
parameters  as 

M3  -  — MS  -  P) 


A1A2  >4i.42 

When  the  shell  deforms,  the  particle  that  had  position  vector  r  in  the  undeformed  state  now 
has  position  vector  R  in  the  deformed  shell.  The  triad  bj  rotates  to  be  Bj.  The  rotation  relating 
these  two  triads  can  be  arbitrarily  large  and  represented  in  the  form  of  a  matrix  of  direction  cosines 
C{xa)  SO  that 

Bj  =  Cijhj  Cij  =  Bj  •  bj  (6) 


A  definition  of  the  2-D  generalized  strain  measures  is  needed  for  the  purpose  of  formulating 
this  problem  in  an  intrinsic  form.  Following  [13]  and  [18],  they  can  be  defined  as 


R.o  =  Aa(Ba  +  -|-  27Q3B3) 


(7) 


and 

Bj^a  =  Ao(  — /Ca2Bl  +  7i'QlB2  +  A’n3B3)  X  Bj  (8) 

where  e^a  are  the  2-D  in-plane  strains,  and  F<ij  are  the  curvatures  of  the  deformed  surface,  which 
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are  the  summation  of  curvatures  of  undeformed  geometry  kij  and  curvatures  introduced  by  the 
deformation  Kij,  and  7^3  are  the  transverse  strains  because  B3  is  not  normal  to  the  reference  surface 
after  deformation.  Please  note  that  the  2-D  generalized  strain  measures  are  defined  by  Eqs.  (7)  and 
(8)  in  an  intrinsic  fashion,  the  symmetry  of  the  inplane  strain  measures  such  that  612  =  £21  does 
not  hold  automatically.  Nevertheless,  one  is  free  to  set  ei2  =  621.  i  e. 

•  ^,2  ^  B2  •  R,l 

A2  ^  ^ 

which  is  a  constraint  used  in  [17]  to  make  the  3-D  formulation  unique. 

At  this  point  sufficient  preliminary  information  has  been  obtained  to  develop  a  geometrically 
nonlinear  shell  theory. 


Compatibility  Equations 

It  is  well  known  that  a  rigid  body  in  3-D  space  has  only  six  degrees  of  freedom.  Thus,  the 
kinematics  of  an  element  of  the  deformed  shell  reference  surface  can  be  expressed  in  terms  of  at 
most  six  independent  quantities:  three  measures  of  displacement,  say  u  •  bj,  and  three  measures  of 
the  rotation  of  B,  (since  the  global  rotation  tensor  C,  which  brings  bf  into  Bj,  can  be  expressed 
in  terms  of  three  independent  quantities).  However,  we  have  the  eleven  2-D  strain  measures  eu, 
2ei2,  €22,  27q3,  Ka3,  and  Kq3  as  defined  in  Eqs.  (7)  and  (8).  Thus,  they  are  not  independent;  there 
are  some  compatibility  equations  among  these  eleven  quantities.  In  [19]  and  [13]  appropriate 
compatibility  equations  are  derived  by  first  enforcing  the  equalities 

R 12  =  R  21  (10) 


and 

Bi,i2  =  B,,21  (11) 

These  two  vector  equations  lead  to  six  independent  compatibility  equations  equivalent  to  a  fomt 
of  those  found  in  [13].  These  equations  are  rewritten  here  for  convenience  in  the  present  notation. 
First,  from  the  B3  components  of  Eq.  (10),  we  obtain 


/I  ,  \  /I  ,  \  (^22723),!  (-412713), 2  ,  ,  T.'  \ 

(1  -t-  e22)«12  ~  (1  +  eil)«21  —  — ~nr - T^t -  ^12(-f'-22  —  All) 


(12) 


A1A2  A1A2 

Next,  from  the  B^  components  ofEq.  (10)  we  obtain  two  equations  for  a=l  and  2,  respectively,  as 

(A2ei2),i  {^1(1 +  eii)],2 


(1  +  622)  Ai3  —  €i2A'23  — 


4iA2 


AIA2 


—  27i3«2i  -f  2723 A'li 


5  OF  23 


(13) 


(1  +  Cll)-1^(^23  “  C12A13  = 


[.4-2(1  +  g22)].i 
.4i  .42 


('4iei2).2 

■4iA2 


—  27i3A'22  +  2723«12 


Finally,  from  the  three  components  of  Eq.  (1 1)  we  have  nine  identities.  However,  there  are  only 
three  independent  equations,  given  by 


{AiKnh 

A1A2 
(.4i/Ci2),2 
A1A2 
(.42Jv23).1 
.4 1^42 


.(  ^  K13K22  -  K12K2Z  =0 

A\A2 

-H  K^Kn  -  =0 

A1A2 


(■4i/vi3),2 

A1A2 


+  K11K22  ~  ^12^21  — 0 


(14) 


There  are  now  eleven  quantities  which  are  related  by  six  compatibility  equations.  This  means  that 
these  strain  measures  can  be  determined  in  terms  of  otily  five  independent  quantities  -  not  six. 

In  the  process  of  dimensional  reduction  of  [17]  to  find  an  accurate  constitutive  model  for  com¬ 
posite  shells,  the  authors  encountered  the  question  whether  one  should  include  K21  and  K12  as  two 
different  generalized  strain  measures.  This  was  determined  by  the  following  argument.  Let  us 
denote  a  new  twist  measure  2uj  =  K12  +  k‘2i-  From  Eq.  (12)  the  difference  between  K21  and  Kn  can 
be  obtained  as 

«12  -  K21  _  ei2(A22  “  Kn)  +  a;(eii  -  €22) 

2  “  (2  +  en+e22) 

This  difference  is  clearly  0(|t)  or  0{^)  disregarding  the  nonlinear  terms  (s  is  the  order  of  gen¬ 
eralized  strains.  h  is  the  thickness  of  the  shell,  I  is  the  wavelength  of  in-plane  deformation  and  R 
is  the  characteristic  radius  of  shell).  One  can  show  that  it  contributes  terms  that  are  0(^^)  or 
O(^)  to  the  strains.  Clearly,  such  terms  will  not  be  counted  in  a  physically  linear  theory  with 
only  correction  up  to  the  order  of  hf  R  and  {h/iy. 

Eqs.  (13)  can  be  solved  for  the  in-plane  curvatures  K13  and  «23,  and  Eq.  (15)  can  be  used 
to  express  K12  and  K21  in  terms  of  a;.  Now,  using  these  expressions,  one  can  rewrite  the  three 
Eqs.  (14)  entirely  in  terms  of  the  eight  strain  measures  en,  2ei2,  €22,  2713,  2723,  Kn,  2uj,  and 
This  confirms  that  only  five  independent  measures  of  displacement  and  rotation  are  necessary  to 
define  these  strain  measures  as  we  will  demonstrate  conclusively  below  by  deriving  such  measures. 


Global  Displacement  and  Rotation  Variables 

There  is  no  unique  choice  for  the  global  deformation  variables.  For  this  reason,  the  importance 
(not  to  mention  the  beauty)  of  an  intrinsic  formulation  is  widely  appreciated.  On  the  other  hand,  for 
the  purpose  of  understanding  the  displacement  field  more  fully,  for  practical  computational  algo- 
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rithms,  and  for  easy  derivation  of  virtual  strain-displacement  relations,  it  is  expedient  to  introduce 
a  suitable  set  of  displacement  measures. 

The  displacement  measures  we  choose  are  derived  by  expressing  R  in  terms  of  r  plus  a  dis¬ 
placement  vector  so  that 

R(a-i,a:2)  =  r(^’i,a:2)  +  Uibj  (16) 

Differentiating  both  sides  of  Eq.  (16)  with  respect  to  Xa,  and  making  use  of  Eq.  (7),  one  can  obtain 
the  identity 

-f  -I-  27a3B3  =  Uijabj  +  Uika  X  b,  (17) 

where  (  );(»  =  The  above  formula  allows  the  determination  of  the  strain  measures  6^$  and 

'2~,a3  in  terms  of  C,  Ui  and  the  derivatives  of  u,-.  Introducing  column  matrices  it  =  [iti  U2  113]^, 
Cl  =  [1  0  OJ^,  €-2  =[0  1  OJ^,  7i  =  [eii  ei2  2713]^,  and  72  =  [e2i  £22  2o'23j^,  we  can  obtain  the 
following  identity  in  matrix  form 


+  7a  —  C{ea  +  -|-  kail') 


(18) 


where  C  is  the  matrix  of  direction  cosines  from  Eq.  (6)  and  ka  is  defined  in  Eq.  (4). 

Rodrigues  parameters  [20]  can  be  used  as  rotation  measures  to  allow  a  compact  expression  of 
C.  These  are  derived  based  on  Euler’s  theorem,  which  shows  that  any  rotation  can  be  represented 
as  a  rotation  of  magnitude  ©  about  a  line  parallel  to  a  unit  vector  e.  Defining  the  Rodrigues 
parameters  /?,  =  2e  •  bj  tan  (f )  and  arranging  these  in  a  column  matrix  p  =  [pi  p2  psj^,  the 
matrix  C  can  simply  be  written  as 


'  1  +  ^ 


(19) 


where  ( )^j  =  —eijk{  )k-  Let  us  also  denote  the  direction  cosines  of  B3  by 

CM=S3i  +  9i  (20) 

Hodges  [21]  has  shown  that,  given  the  third  row  of  C,  the  Rodrigues  parameters  can  be  uniquely 


7  OF  23 


expressed  in  terms  of  di  as 


Pi 

P2 


pzd\  —  2^2 

2  +  03 
P302  +  201 
2  +  03 


P3  =2  tan 


(21) 


where  ps  can  be  understood  as  a  change  of  variables  to  simplify  later  parts  of  the  derivation.  Later 
on  we  will  discuss  the  meaning  of  <t>z  for  a  special  case.  Finally,  it  is  noted  that  the  three  rotational 
parameters  0;  are  not  independent  but  instead  satisfy  the  constraint 

01  +  02  +  (1  +  03)^  —  1  (22) 


When  Eq.  (21)  is  substituted  into  Eq.  (19),  the  resulting  elements  of  C  can  be  expressed  as 
functions  of  0,  and  03 


On 

C12 

Ci3 

^21 

C22 

<^23 

<^31 

Cz2 

Czz 


(2  +  03  -  0j)  cos dz  -  0102  sino3 


2  +  ^3 

(2  +  03- 

02  )  sin  03  —  0102  cos  G>3 

2  +  ^3 

=  -  01  cos  03  —  02  sin  03 

—  (2  +  03 

-  01 )  sin  03  —  0102  cos  03 

2  +  03 

_  (2  +  03  — 

01)  cos  03  +  0102  sin  03 

2  +  03 

=01  sin  (j)z  —  02  cos  03 

=01 

=02 

=1+03 


(23) 


This  representation  reduces  to  those  of  [22]  when  considering  small,  finite  rotations.  There  is  an 
apparent  singularity  in  the  present  scheme  when  03  =  —2  (Le.,  when  the  shell  deforms  in  such  a 
way  that  B3  is  pointed  in  the  opposite  direction  of  b3).  This  should  pose  no  practical  problem, 
however,  since  0i  =  02  =  0  for  that  condition,  and  none  of  the  kinematical  relations  become 
infinite  in  the  limit  as  03  — >  —2. 

When  these  expressions  for  the  direction  cosines  are  substituted  into  Eq.  (18),  explicit  expres- 
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sions  for  the  strain  measures  can  be  found  as 

=  [(2  +  %  -  ^;)(1  +  “1:1  -  +  feliua)  -  gl»»(“2:l  +  fcl3.».l)  ^ 

L  2  +  03  J 

+  [(2  +  ^3  -  BI){u2,x  +  fciatti)  -  ^1^2(1  +  txi;i  -  fci3tt2  +  ^  02(fcii«i  -  «3;i)l  sin  <^3  -  1 

[  2  +  03  ’  J 

e,.,  =  [(2  +  ^3  -  e^^)il  +  U2;2  +  hzui  +  fc22U3)  -  glg2(ui;2  -  fc23U2)  _ 

■■  [  2  +  ^3  ’  J 

+  [(2  +  ^3  -  0l){h3U2  -  Uia)  +  ^1^2(1  +  U2:2  +  fe23»l  +  ^2^)  ^  ^  ^ 

[  2  +  03  ’  J 

,,,  =  (2  +  .^3  -  g|)(»2:l  +  fcl3»l)  -  ^1^2(1  +  ni;i  -  fci3»2  +  knUz)  J 

2  +  03  '  J 


(2  +  ^3  -  gf)(l  +  uui  -  knU2  +  fciiU3)  -  O102{U2-1  +  fci3_»0  +  0^(knui  -  U3.i)l  sin<?3 

.  2  +  03  ■  J 


_  (2  +  03  —  0i)(»l:2  ~  fc23»2)  +  01^2(1  +  U2:2  +  fcasUl  +  ^22^3)  2  —  ^22^2) 

2  +  03 

+  (2  +  03  -  0i)(l  +  U2:2  +  hsui  +  k22U3)  -  0lg2(ui:2  -  ^23^2)  _  ^  -  k22U2)  \  sin  63 

2  +  03 


cos  03 


2^.13  =  0i(l  +  14;1  —  A"i3U2  +  fciiU3)  +  02C^2;1  +  +  (1  +  ('^^3:1  ”  knlti) 

”  A’23W2)  +  ^2(1  +  ^2:2  +  +  ^22^3)  +  (1  +  ^s)  ('^*3:2  ~  ^^22^^2) 


These  expressions  explicitly  depend  on  sin  ^>3  and  cos  <63.  It  is  evident  that  one  can  choose  63  so 
that  612  =  621,  yielding 

_  ni  +  0§(»2:1  +  ki3Ui)  -  0i(ui;2  -  k23U2)  +  0102  [«1:1  -  »2:2  +  (^11  "  ^22)^3  -  ^23^1  -  A-'isH 

tan  03  -  ^  02(uj.j  4-  knU3  -  fcl3’“2)  +  ^2(^2;2  +  h'l’-i'S  +  hiU^l)  +  ^1^2(ui;2  +  U2;l  “  ^‘23f‘2  +  kuUi) 

(25) 

where 


ni  =  (2  +  03)  [ui;2  -  •U2;l  -  0l(U3;2  -  k22'tJ'2)  +  02(ti3;l  “  ^H'Wi)  “  hiUi  -  ^23 “2] 
n2  =  (2  +  03)[0l(W3;l  -  hl^l)  +  02{u3-2  -  ^22^2)  -  Ml;l  “  ■l^2;2  —  2  —  63 
—  k23Ui  —  {k22  +  A:ii)u3  +  /C13W2] 

It  is  now  clear  that  once  the  functions  ui,  U2,  U3,  0i  and  62  are  known,  the  entire  deformation  is 
determined.  Because  of  this,  one  should  expect  that  a  variational  formulation  would  yield  only  five 
equilibrium  equations  -  no/ 3'/.v. 

For  small  displacement  and  small  strain,  one  can  obtain  ^3  as 

_  (■^2/*2),i  (■^rt^i),2  (27) 

~  2A1A2 
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which  is  the  angle  of  rotation  about  B3,  the  same  as  obtained  in  [23]. 

Although  one  can  now  find  exact  expressions  for  eu,  2ei2,  €22.  2713  and  2723  which  are  inde¬ 
pendent  of  ^3,  such  expressions  are  rather  lengthy  and  are  not  given  here.  Alternatively,  one  could 
leave  ^3  in  the  equations  and  regard  Eq.  (25)  as  a  constraint.  This  would  allow  the  construction  of 
a  shell  finite  element  which  would  be  compatible  with  beam  elements  which  have  three  rotational 
degrees  of  freedom  at  the  nodes. 

Expressions  for  the  curvatures  can  be  found  in  terms  of  C  as 

IQ  =  -C.jf'  +  (28) 


where 

Ha  =  [— ^'al  «al 

Following  [24],  the  curvature  vector  can  also  be  found  using  Rodrigues  parameters 

I  -I 

Ha  =  - - ’—^P.a  +  Cka 

1  A.  £i£ 


(29) 


(30) 


Using  the  form  of  C  from  Eqs.  (23),  the  curvatures  become 


where 


^  .  ^3:o(^lCOS03  +  ^2Smd'3)  ,  ,  , 

Kal  =  “ha  COS  C>3  +  02:q  Sill  <?3 - ^rTTi - b  Kal  —  KqI 


Kq2  =  ~^l-,a  sin  03  +  ^2;q  COS  03  -f 
&ha&2  —  ^l&2\a 


2  -F  ^3 

^3:a(^l  sin  03  -  $2  COS  03) 

2  +  ^3 


+  ka2  —  ka2 


K^aS  —  C>3:a  + 


2  -f-  ^3 


+  k 


'^a3 


(31) 


kal  —  (^a3 
ka2  (^0:3 


ka2^l  kf^\02 


+  o  T  J')(^i  sin  03  -  O2  cos  03)  -1-  kc2  sin  03  +  kd  cos  03 


2  -F  ^3  2  -F  ^3 

ka2^l  ,  kd02 


2  4-  ^3  2  -F  ^3 

ka3  =  ~ka2^1  +  kal02  +  ^a3^3 


+  "  ^  )(^1  cos  03  +  02  sin  03)  +  kal  COS  03  -  kd  sill  ©3 


(32) 


As  before,  03  can  be  eliminated  from  these  expressions,  so  that  all  six  curvatures  can  be  expressed 
in  terms  of  five  independent  quantities.  Note  that  Kd  are  not  independent  2-D  generalized  strains. 
They  will,  however,  appear  in  the  equilibrium  equations  because  of  their  appearance  in  the  virtual 
strain-displacement  relations  derived  later. 
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2-D  Constitutive  Law 

To  complete  the  analysis  for  an  elastic  shell,  a  2-D  constitutive  law  is  required  to  relate  2- 
D  generalized  stresses  and  strains.  As  mentioned  before  the  constitutive  law  can  not  be  exact, 
however,  one  should  try  to  avoid  introducing  any  unnecessary  approximation  in  addition  to  the 
already-approximate  3-D  constitutive  relations. 

Among  many  approaches  that  have  been  proposed  to  deal  with  dimensional  reduction,  the  ap¬ 
proach  in  [17]  stands  out  for  its  accuracy  and  simplicity.  In  that  work,  a  simple  Reissner-Mindlin 
type  energy  model  is  constructed  that  is  as  close  as  possible  to  being  asymptotically  correct.  More¬ 
over,  the  original  3-D  results  can  be  recovered  accurately.  The  resulting  model  can  be  expressed 
as 

2n  =  €^Ae+n^G')+2€^F  (33) 

where  e  =  [cn  2ei2  622  kh  «i2  +  ><21  «22j^  and  7  =  [2713  2723]^.  It  is  noticed  that  there  is 
only  one  in-plane  shear  strain  612  in  Eq.  (33).  This  is  possible  only  after  one  uses  the  constraints 
in  Eq.  (9).  Moreover,  the  strain  energy  is  independent  of  so  that  the  rotation  about  the  normal 
only  appears  algebraically,  making  it  possible  for  it  to  be  eliminated. 

This  simple  constitutive  model  is  rigorously  reduced  from  the  original  3-D  model  for  multilayer 
shells,  each  layer  of  which  is  made  with  an  anisotropic  material  with  monoclinic  symmetry.  The 
variational  asymptotic  method  [16]  is  used  to  guarantee  the  resulting  2-D  shell  model  to  yield  the 
best  approximation  to  the  energy  stored  in  the  original  3-D  structure  by  discarding  all  the  insignif¬ 
icant  contribution  to  the  energy  higher  than  the  order  of  {h/lY  and  h/R.  The  stiffness  matrices 
.4  and  G  obtained  through  this  process  carry  all  the  material  and  geometry  information  through 
the  thickness  (see  Eqs.  (63)  and  (73)  in  Ref.  [17]  for  detailed  expressions).  The  term  containing 
the  column  matrix  F  is  produced  by  body  forces  in  the  shell  structure  and  tractions  on  the  top  and 
bottom  surfaces,  and  it  is  very  important  for  the  recovery  of  the  original  3-D  results.  Interested 
readers  can  refer  to  Ref.  [17]  for  details  of  constructing  the  model  in  Eq.  (33)  for  multilayered 
composite  shells. 

Having  obtained  the  2-D  constitutive  law  from  3-D  elasticity,  one  can  derive  all  the  other 
relations  oyer  the  reference  surface  of  the  shell,  a  2-D  continuum. 

Virtual  Strain  Displacement  Relations 

In  order  to  derive  intrinsic  equilibrium  equations  from  the  2-D  energy,  it  is  necessary  to  express 
the  variations  of  generalized  strain  measures  in  terms  of  virtual  displacements  and  virtual  rotations. 
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The  variation  of  the  energy  expressed  in  Eq.  (33)  can  be  written  as 


d'n  = 


+ 


an  . 

<)eii  + 


an 


aeii 

an 


d' 


m 


^723  + 


aei2 

an 


5ei2  + 


an , 

^€22  + 


an 


ae22 


57: 


dm 

,  an^  an  . 

7^ - 5k11  +  T^ou  +  7: - OK22 

a^ii  au  OK22 
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(34) 


It  is  now  obvious  that  one  must  express  5eii,  . . SK22,  in  terms  of  virtual  displacements  and 
rotations  in  order  to  obtain  the  final  Euler-Lagrange  equations  of  the  energy  functional  in  their 
intrinsic  form.  Following  Ref.  [24],  we  introduce  measures  of  virtual  displacement  and  rotation 
that  are  “compatible”  with  the  intrinsic  strain  measures.  For  the  virtual  displacement,  we  note  the 
form  of  Eq.  ( 1 8)  and  choose 

Jq  =  C5u  (35) 


Similarly,  for  the  virtual  rotation,  we  note  the  form  of  Eq.  (28)  and  write 


5i)  = 


(36) 


where  6v  is  a  column  matrix  arranged  similarly  as  the  curvature  column  matrix  in  Eq.  (4)  6v  = 
[— 5r'.,  ^3]^.  The  bars  indicate  that  these  quantities  are  not  necessarily  the  variations  of 

functions.  Using  these  relations  it  is  clear  that 

5u  =  (37) 


and 

SC  =  -JuC  (38) 


Let  us  begin  with  the  generalized  strain-displacement  relationship,  Eq.  (18).  A  particular  in-plane 
strain  element  can  be  written  as 


=  ej  C{ea  +  U;a  +  AJqU)  —  Cq 


(39) 


Taking  a  straightforward  variation,  one  obtains 

=  ej  [^^(eQ  -f-  U.^a  +  kaU)  +  C{6iI-ol  +  ka^u) 


(40) 


The  right  hand  side  contains  and  Su-a^  which  must  be  eliminated  in  order  to  obtain  variations  of 
the  strain  that  are  independent  of  displacements.  These  are  needed  to  derive  intrinsic  equilibrium 
equations. 

Premultiplying  both  sides  of  Eq.  (18)  by  C'^,  making  use  of  Eq.  (36),  and  finally  using  a 
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property  of  the  tilde  operator  that,  for  arbitrary  column  matrices  V  and  Z,  YZ  =  —ZY,  one  can 
make  the  first  term  in  brackets  on  the  righthand  side  independent  of  U;q.  After  all  this,  one  obtains 


5C{ea  +  u-a  +  kaU)  =  5CC^  {e^  +  la)  =  -H{ea  +  7a)  =  {^a  +  la)^'^  (41) 


An  expression  for  the  second  term  in  brackets  on  the  right  hand  side  of  Eq.  (40)  can  now  be 
obtained  by  differentiating  Eq.  (37)  with  respect  to  Xq  and  premultiplying  by  C.  This  yields 


C{5u.oc  +  kaSu)  —  C{C^6q);a  +  CkocSu  =  Sq.^  +  KaSq  (42) 

Substituting  Eqs.  (41)  and  (42)  into  Eq.  (40),  one  obtains  an  intrinsic  expression  for  the  variation 
of  the  in-plane  strain  components  as 


d€a3  = 


^Q:a  +  +  (Ca  + 


(43) 


where  e^e'a  vanishes  when  a  =  This  matrix  equation  can  be  written  explicitly  as  four  scalar 
equations 


^en  =  —  Ai3^?2  +  ~  27i3cib'i  -f 

Sei2  —  <5g2:i  +  KizSqi  -i-  ki2^93  ~  -ln^^'2  ~  (1-  + 

^€21  =  5Qi.2  —  K2z5(l2  +  /^2I^93  ~  “^IZZ^Wl  +  (1  +  ^22)^^’3 

5e22  =  ^q2;2  +  +  ^^•22^Q3  ~  ‘^l2Z^4-'2  ~  ^12(^^3  (44) 


The  variations  6ei2  and  ^€21  should  be  equal  due  to  Eq.  (9);  hence,  one  can  solve  for  the  virtual 
rotation  component  about  B3  as 


<^921  ~  <^9i-2  “I"  +  ^23^0.2  +  («:i2  ~  ^2l)^<lz  ~  2713^02  +  ‘^l2Z^Wl  . 

=  — : - - - - - -  (45) 

Z  +  eii  +  €22 


It  is  now  possible  to  write  the  variations  of  all  strain  measures  in  terms  of  three  virtual  displacement 
and  two  virtual  rotation  components  as 


5eii  =  —  K\i6q2  +  KnSq^  —  2iizSipi  -h  eizStp^ 

S€22  =  ^92:2  +  ^'^2zSqi  +  /V22d'93  "  ‘^IzzH'z  ~  ^nH'z 

2^612  =  5g2;i  +  ^9i;2  +  ~  I^zz^Qi  +  2u;d'93 

~  27i35'^2  “  2723<^'|/’i  +  (^22  ~  fll)<^'03  (46) 

with  taken  from  Eq.  (45). 
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Let  us  now  consider  the  transverse  shear  strains 


-  ^3  tljQ:  ATofU)  6q; 


(47) 


Following  a  procedure  similar  to  the  above,  one  can  obtain  the  virtual  strain-displacement 
equation  for  transverse  shear  strains  as 


257q3  =  63  -f  Ka^q  +  (ea  +  7q)<^0 


(48) 


Explicit  expressions  for  the  variations  of  the  shear  strain  components  are  now  easily  written  as 


257q3  =  "b  ^(x3^W^  ~ 


(49) 


Filially,  variations  of  the  curvatures  are  found.  First,  taking  the  straightforward  variation  of  Eq.  (28), 
one  obtains 

S7C  =  -  +  dCZc^  +  cTjC'^  (50) 


A. 


Ac 


In  order  to  eliminate  <5C.q,  we  differentiate  Eq.  (36)  with  respect  to  Xa 


=  -SCaC^  -  SCCl 


(51) 


In  order  to  eliminate  SC,  we  can  use  Eq.  (38).  Then,  Eq.  (50)  becomes 


SKc  =  Slp.^  +  KaSp  -  SlpKa 


(52) 


Using  another  tilde  identity  (V'Z  =  V' Z  -  ZY)  one  can  find  the  virtual  strain-displacement  relation 
as 

SKcc  =  +  KaSlp  (53) 

In  explicit  form 


Sipii  — 

<5aCii  —  — - Y 

Ai 

S'lpo  O  '  " 

Sk>22  =  — - 1“  ”  ^21<^03 

A2 

25uj  =  -I-  +  A"i3^i  -  /V23b'0.2  +  {K'n  - 

A2  .4i 


(54) 


where  Sii’n  can  again  be  elirhinated  by  using  Eq.  (45). 
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Intrinsic  Equilibrium  Equations 

In  this  section,  we  will  make  use  of  the  virtual  strain-displacement  relations  in  the  variation  of 
the  internal  strain  energy  in  order  to  derive  the  intrinsic  equilibrium  equations.  Here  we  define  the 
generalized  forces  as 


To  use  the  principle  of  virtual  work  to  derive  the  equilibrium  equations,  one  needs  to  know 
the  applied  loads.  In  addition  to  the  applied  loads  used  in  the  modeling  process,  r,Bi  at  the  top 
surface,  ,JjBi  at  the  bottom  surface  and  body  force  OiB,-  [17],  one  can  also  specify  appropriate 
combinations  of  displacements,  rotations  (geometrical  boundary  conditions),  running  forces  and 
moments  (natural  boundary  conditions)  along  the  boundary  around  the  reference  surface.  It  is 
trivial  to  apply  the  geometrical  boundary  conditions.  Although  it  is  possible  in  most  cases  that 
natural  boundary  conditions  can  be  derived  from  Newton’s  law,  the  procedure  is  tedious  and  not 
easily  applied  here  because  the  physical  meanings  for  some  of  the  generalized  forces  are  not  clear. 
Thus,  natural  boundary  conditions  are  best  derived  from  the  principle  of  virtual  work. 

Suppose  on  boundary  F  (see  Fig.  2),  we  specify  a  force  resultant  Nuu  and  moment  resultant 
Muu  along  the  outward  normal  of  the  boundary  curve  tangent  to  the  reference  surface  v,  Ni,t  and 
along  the  tangent  of  the  boundary  curve  r,  along  the  normal  of  the  reference  surface. 
Then  the  principle  of  virtual  work  (strictly  speaking,  the  principle  of  virtual  displacements)  can  be 
stated  as: 


J  J{5U  -  5qji  -  Stpa'>^a)AiA2dxidx2- 


I  Nut^Qt  +  A/i/r^'07-)cir  0 


where  /j  and  ma  are  taken  directly  from  [17]. 

It  is  now  possible  to  obtain  intrinsic  equilibrium  equations  and  consistent  edge  conditions  by 
use  of  the  principle  of  virtual  work  and  the  virtual  strain-displacement  relations  derived  in  the 
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previous  section.  The  equilibrium  equations  are 

(A2iVii).i  ,  [A,{Nn+^f)\2 


—  Ki3{Ni2  —  M)  —  K23N22  +  Qlt^n  +  Q2^21  +  /l  =  0 
+  K23{Ni2  +  a/")  +  +  QiK.12  +  Q2^22  +  /2  =  0 


A1A2  A1A2 

iAiN22),2  ^  [A2iNx2-Af)\, 

-4i-42  A1A2 

i — j*— ^  ^  ~  ^"^11 ‘^11  K22N22  ““  2a.’A^i2  +  (/^i2  “  ^2\)Af  +  /s  =  0 

-■11.42  ^^1^2 

(.4,.U„),  ,  (A,Mn)x  ^  , 

—  Q2fi2  +  27i3iVii  +  2723(^^12  +  A/*)  —  AI12K13  —  M22K23  +  nil  —  0 

(A24^/i2).1  ,  (Aiil/22),2 


+ 


—  Q2(1  +  £22) 


A1A2  .■41^2 

—  <5x^12  +  27i3(iVi2  ~  A/")  +  27'23A'22  +  i\/ii ATn  +  M12K23  +  rni  —  0  (57) 


where 

_  {N22  —  A'ii)6i2  +  Ni2{€n  —  €22)  +  Al22f^21  ~  ■‘^Al^l2  +  A/i2(A"ll  —  A22) 

2  +  eii  +  e22 

The  associated  natural  boundary  conditions  on  F  are 

Ni/u  =  +  2nin2N\2  +  'it^N22 

N^t  =  nin2{N22  ~  A^ii)  +  (^i  ~  n\)Ni2  —  M 
Nt/u  ~  A^ii  +  2nin2Ni2  *b  v^N22 
Av3  =  ttiQi  +  n2Q2 
AIuv  —  tijA/ji  +  2nin2Mi2  •\-‘tT^M22 
Mut  =  nin2{M22  ~  A/ll)  +  (t^i  ~  ti2)A/i2 


(58) 


(59) 


where  rii  =  cos  0,  n2  =  sin  0,  and  (f)  is  the  angle  between  the  outward  normal  of  the  bbundary 
and  xi  direction  as  shown  in  Fig.  2.  The  terms  containing  A/  stem  from  consistent  inclusion  of  the 
finite  rotation  from  undeformed  triad  to  deformed  triad  although  the  nonzero  rotation  about  B3  is 
expressed  in  terms  of  other  kinematical  quantities.  Similar  terms  are  found  in  the  shell  equations 
derived  by  Berdichevsky  [16]  where  only  five  equilibrium  equations  are  derived. 

In  a  mixed  formulation,  A/*  can  be  shown  to  be  the  Lagrange  multiplier  that  enforces  Eq.  (45). 
To  further  understand  the  nature  of  A/”  one  can  undertake  the  following  exercise:  Setting  Pi  =  0 
and  ei2  =  £21  for  the  equilibrium  equations  given  in  [13],  can  be  solved  from  Reissner’s 

sixth  equilibrium  equation.  This  shows  that  Reissner’s  is  the  same  as  our  Af,  and  Reiss- 
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net’s  is  the  same  as  our  N^-  Finally,  substitution  of  this  sixth  equation  into  the  other  five 

yields  the  five  equilibrium  equations  given  here  in  Eqs.  (57).  It  is  noted  that  Reissner’s  equilib¬ 
rium  equations  are  derived  based  on  the  basis  of  Newton’s  law  of  motion  without  consideration  of 
either  constitutive  law  or  strain-displacement  relations.  However,  the  present  derivation  is  purely 
displacement-based.  The  reproduction  of  those  equilibrium  equations  by  the  present  derivation 
illustrates  that,  as  long  as  the  formulation  is  geometrically  exact,  one  can  derive  exact  equilibrium 
equations. 

A  few  investigators  have  noted  an  apparent  conflict  between  the  symmetry  of  the  stress  resul¬ 
tants  and  the  satisfaction  of  moment  equilibrium  about  the  normal.  In  reality  there  is  no  conflict, 
but  one  must  be  careful.  We  have  shown  herein  that  the  triad  B,  can  always  be  chosen  so  that 
ei-2  =  If  this  relation  is  enforced  strongly,  there  is  only  one  in-plane  shear  stress  resultant, 
:Vi2.  that  can  be  derived  from  the  energy.  In  that  case  the  physical  quantity  associated  with  the 
antisymmetric  part  of  Reissner’s  in-plane  stress  resultants,  while  it  is  not  available  from  the  con¬ 
stitutive  law,  is  nevertheless  available  as  a  reactive  quantity  from  the  moment  equilibrium  equation 
about  the  normal.  However,  it  must  be  stressed  that  the  moment  equilibrium  equation  about  the 
normal  is  not  available  from  a  conventional  energy  approach,  in  which  the  virtual  displacements 
and  rotations  must  be  independent. 

In  a  somewhat  similar  vein,  not  being  able  to  obtain  the  antisymmetric  part  of  the  moment  stress 
resultants  from  derivatives  of  the  2-D  strain  energy  is  a  result  of  the  approximate  dimensional  re¬ 
duction  process  in  which  it  was  determined,  based  on  asymptotic  considerations  and  geometrically 
nonlinear  3-D  elasticity,  that  the  antisymmetric  term  —  K21  does  not  appear  as  an  independent 
generalized  strain  measure  in  the  2-D  constitutive  law  with  correction  only  to  the  order  of  h/R. 
However,  if  a  more  refined  theory  with  respect  to  h/R  is  required,  then  K12  —  K21  would  appear  as 
a  generalized  strain  in  the  2-D  constitutive  law  and  a  new  generalized  moment  would  be  defined 
based  on  the  constitutive  law. 

For  practical  computational  schemes,  equilibrium  equations  and  boundary  conditions  need  to 
use  the  constitutive  law  to  relate  with  the  generalized  2-D  strains.  Finally  a  set  of  kinematical 
equations  is  needed.  Depending  on  how  this  part  is  done,  the  analysis  can  be  Completed  in  either 
of  two  fundamentally  different  ways:  a  purely  intrinsic  form,  relying  on  compatibility  equations, 
and  a  mixed  form  relying  on  explicit  strain-displacement  relations. 

In  the  intrinsic  form  we  have  five  equilibrium  equations,  Eqs.  (57);  six  compatibility  equations, 
Eqs.  (12)  — (14);  and  the  eight  constitutive  equations  -  a  total  of  19  equations.  The  19  unknowns  are 
the  eight  stress  resultants,  NnyNn,  N22,  Qi,  Q2,  Mi.  Mn,  and  M22;  and  the  1 1  strain  measures 
Cii,  2612,  ^22.  2''/i3,  2')'’23.  Atii.  2ci2i2>  and  ^22.  along  with  kj3,  /c23>  and  K12  ^21  •  The  last  three  strain 

rheasures  appear  in  the  equilibrium  equations  but  not  in  the  constitutive  law. 

In  a  mixed  formulation  one  would  use  the  same  five  equilibrium  equations  and  eight  con- 


17  OF  23 


stitutive  equations.  One  would  also  need  a  set  of  strain-displacement  relations  among  the  1 1 
generalized  strain  measures  cn,  2ei2,  £22.  2713,  2723,  /tii.  2u;,  and  K22,  along  with  K13,  K23,  and 
kv2  -  «2i,  and  the  five  global  displacement  and  rotational  variables  ui,  U2,  U3,  61,  and  62.  One 
possible  set  of  such  equations  is  as  follows:  use  five  of  Eqs.  (24),  using  either  €12  or  e2i;  use  the 
six  Eqs.  (31).  There  are  also  the  two  other  rotational  variables  63  and  (p3,  which  are  governed 
by  Eqs.  (22)  and  (25),  respectively.  This  way  there  are  26  equations  and  26  unknowns.  This 
mixed  formulation  is  capable  of  handling  boundary  conditions  on  2-D  stress  resultants  and  dis¬ 
placement/rotation  variables.  At  least  in  principle,  one  could  recover  a  displacement  formulation 
by  eliminating  all  the  unknowns  except  the  displacement  and  rotation  variables. 

Eqs.  (57)  and  (58)  contain  terms  that  could  be  disregarded  because  of  the  original  assumption 
of  small  strain.  We  will  not  undertake  this  simplification  here,  because  it  is  out  of  the  scope  of  the 
present  study  to  actually  implement  the  2-D  nonlinear  theory.  Therefore,  our  equilibrium  equations 
and  kinematical  equations  are  geometrically  exact;  all  approximations  stem  from  the  dimensional 
reduction  process  used  to  obtain  the  2-D  constitutive  law. 

The  present  work  is  a  direct  extension  of  [  18]  to  treat  shells.  If  one  sets  kij  =  0  and  .4^  =  1,  all 
the  formulas  developed  here  will  reduce  to  those  in  [18],  which  indirectly  verifies  that  derivation. 

Conclusions 

A  nonlinear  shear-deformable  shell  theory  has  been  developed  to  be  completely  compatible 
with  the  modeling  process  in  [17].  The  compatibility  equations,  kinematical  relations  and  equi¬ 
librium  equations  are  derived  for  arbitrarily  large  displacements  and  rotations  under  the  restriction 
that  the  strain  must  be  small.  The  resulting  formulae  are  compared  with  others  in  the  literature. 
The  following  conclusions  can  be  drawn  from  the  present  work: 

1 .  The  variational  asymptotic  method  can  be  used  to  decouple  the  original  3-D  elasticity  prob¬ 
lem  of  a  shell  into  a  l-D,  through-the-thickness  analysis  [17]  and  a  2-D,  shell  analysis.  The 
through-the-thickness  analysis  provides  both  an  accurate  2-D  constitutive  law  for  the  nonlin¬ 
ear  shell  theory  and  accurate  through-the-thickness  recovery  relations  for  3-D  displacement, 
strain,  and  stress.  This  way,  an  intimate  relation  between  the  shell  theory  and  3-D  elasticity 
is  established. 

2.  A  full  finite  rotation  must  be  applied  to  fully  specify  the  displacement  field.  However,  since 
the  strain  energy  on  which  the  formulation  is  based  is  independent  of  /?q3,  the  rotation  about 
the  normal  is  not  independent  and  can  be  expressed  in  terms  of  other  quantities.  Thus,  it 
can  be  chosen  so  that  the  two-dimensional,  in-plane  shear  strain  measures  are  equal.  This 
way  all  the  strain  measures  can  be  expressed  in  terms  of  five  independent  quantities:  three 
displacement  and  two  rotation  measures,  and  only  one  stress  resultant  for  in-plane  shear  can 
be  derived  from  the  2-D  energy. 
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3.  Only  five  equilibrium  equations  are  obtainable  in  a  displacement-based  variational  formula¬ 
tion.  Moment  equilibrium  about  the  normal  is  satisfied  implicitly.  If  one  does  not  include 
the  full  finite  rotation,  but  rather  sets  the  rotation  about  the  normal  equal  to  zero,  the  cor¬ 
rect  equilibrium  equations  cannot  be  obtained.  This  should  shed  some  light  on  the  nature  of 
“drilling”  degrees  of  freedom. 
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Appendix  E 


An  Asymptotic  Approach  for 
Thermoelastic  Analysis  of  Laminated 

Composite  Plates 

Wenbin  Yu*  and  Dewey  H.  Hodges^ 

Georgia  Institute  of  Technology,  Atlanta,  Georgia  30332-0150 

A  thermoelastic  model  for  analyzing  laminated  composite  plates  under  both  mechanical 
and  thermal  loadings  is  constructed  by  the  variational  asymptotic  method.  The  original 
three-dimensional  nonlinear  thermoelasticity  problem  is  formulated  based  on  a  set  of  intrinsic 
\'ariables  defined  on  the  reference  plane  and  for  arbitrary  deformation  of  the  normal  line. 
Then  the  variational  asymptotic  method  is  used  to  rigorously  split  the  three-dimensional 
problem  into  two  problems:  a  nonlinear,  two-dimensional,  plate  analysis  over  the  reference 
plane  to  obtain  the  global  deformation  and  a  linear  analysis  through  the  thickness  to  provide 
the  two-dimensional  generalized  constitutive  law  and  the  recovering  relations  to  approximate 
the  original  three-dimensional  results.  The  nonuniciueness  of  asymptotic  theory  correct  up 
to  a  certain  order  is  used  to  cast  the  obtained  asymptotically  correct  second-order  free 
energy  into  a  Reissner-Mindlin  type  model  to  account  for  transverse  shear  deformation.  The 
present  theory  is  implemented  into  the  computer  program,  Variational  Asymptotic  Plate 
and  Shell  Analysis  (VAPAS).  Results  from  VAPAS  for  several  cases  have  been  compared 
with  the  exact  thermoelasticity  solutions,  classical  lamination  theory  and  first-order  shear- 
deformation  theory  to  demonstrate  the  accuracy  and  power  of  the  proposed  theory. 

Keywords:  Asymptotic  series,  anisotropic  plates,  finite  element  method,  strain  distribution, 
stress  distribution,  thermoelasticity 

Introduction 

Composite  materials  have  found  increasing  applications  in  engineering  practices  due  to 
their  superior  engineering  properties  and  improving  manufacturing  technology.  However,  the 
heterogeneity  and  anisotropy  of  such  materials  make  the  traditional  analysis  method  used 
for  designing  homogeneous  and  isotropic  structures  obsolete.  Moreover,  structures  made 
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■with  composite  materials  are  more  sensitive  and  vulnerable  to  temperature  change  than 
their  isotropic  counterpart.  The  reason  is  that  the  thermal  expansion  coefficients  of  different 
constituents  of  the  material  are  usually  dramatically  different  from  each  other  resulting  in 
high  stresses  due  to  sudden  temperature  change.  The  analysis  including  thermal  effects  is 
much  more  involved  than  that  for  isothermal  conditions. 

Many  engineering  structures  made  with  composite  materials  have  one  dimension  much 
smaller  than  the  other  two  and  can  be  modeled  as  plates.  Only  a  few  exact  solutions  exist  for 
very  idealized  cases  (see  Savoia  and  Reddy  (1995)  and  the  references  cited  there).  Researchers 
are  trying  to  develop  simplified  models  to  provide  an  approximate  representation  for  more 
general  cases.  Within  the  last  few  decades,  a  tremendous  research  effort  has  been  invested 
in  this  area,  and  various  approximate  models  have  been  proposed  (Wu  and  Tauchert  (1980); 
Xoor  and  Burton  (1992);  Reddy  (1997);  Noor  and  Malik  (2000);  Rohwer  et  al.  (2001)). 
Generally  speaking,  these  models  are  derived  from  three-dimensional  (3-D)  thermoelasticity 
theory,  making  use  of  the  fact  that  the  plate  is  thin  in  some  sense.  Although  it  is  plausible 
to  consider  the  smallness  of  the  thickness  of  plate  structures,  construction  of  an  accurate 
two-dimensional  (2-D)  model  for  a  3-D  body  still  introduces  a  lot  of  challenges.  Almost  all 
the  proposed  models  in  the  literature  can  give  a  good  prediction  of  the  global  behavior  of 
the  plate.  However,  they  have  serious  difficulties  in  providing  accurate  distributions  of  the 
3-D  displacements,  strains,  and  stresses  through  the  thickness.  Part  of  the  reason  is  that 
most  models  adopt  ad  hoc  assumptions  (such  as  having  displacement  or  stress  components 
vary  through  the  thickness  according  to  a  certain  function)  which  violate  the  exact  solutions. 
For  example,  most  high  order  theories  (except  for  some  layerwise  theories  such  as  Cho  and 
Averill  (2000))  assume  a  continuity  for  the  3-D  displacement  field  through  the  thickness 
which  in  reality  they  are  piecewise  continuous.  In  the  case  of  thermal  loading,  the  prediction 
of  these  ad  hoc  models  become  even  worse,  if  not  wrong,  and  the  results  should  be  examined 
more  cautiously. 

It  must  be  understood  that  all  plate  theories,  no  matter  how  involved  they  may  appear, 
are  inherently  approximate-  The  approximation  lies  in  the  2-D  constitutive  law  relating 
2-D  strains  and  stress  resultants,  which  is  a  direct  consequence  of  eliminating  the  thickness 
coordinate  from  the  independent  variables  of  the  governing  equations  of  the  boundary-vnlue 
problem  for  a  plate.  This  sort  of  approximation  is  inevitable  if  one  wants  to  take  advantage 
of  the  smallness  of  the  thickness  to  simplify  the  analysis.  It  is  interesting  to  note  that  the  3-D 
constitutive  relations  are  essentially  approximate  and  determined  by  experiments.  However, 
this  cannot  be  used  as  an  excuse  to  introduce  unnecessary  assumptions.  For  example,  for 
small-strain  analysis  of  plates,  it  is  reasonable  to  assume  that  the  thickness,  h,  is  small 
compared  to  the  wavelength  of  deformation  of  the  reference  surface,  1.  However,  it  is  not  at 
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all  reasonable  to  assume  a  priori  some  ad  hoc  displacement  field,  although  that  is  the  way 
most  plate  theories  are  constructed. 

In  this  paper,  we  first  cast  the  original  3-D  thermoelasticity  problem  in  an  intrinsic 
form  so  that  the  theory  is  applicable  for  arbitrarily  large  displacement  and  global  rotation, 
subject  only  to  the  strain’s  being  small;  see  Hodges  et  al.  (1993).  Then,  the  Variational 
Asymptotic  Method  (VAM),  introduced  by  Berdichevsky  (1979),  is  used  to  split  the  original 
nonlinear  3-D  elasticity  problem  into  a  linear,  one-dimensional  (1-D),  through-the-thickness 
analysis  and  a  nonlinear,  2-D,  plate  analysis.  The  through-the-thickness  analysis  produces 
the  2-D  constitutive  law  and  non-mechanical  stress  resultants  to  be  used  in  the  2-D  plate 
analysis,  along  with  recovering  relations  that  yield  the  3-D  displacement,  strain  and  stress 
fields  through  the  thickness  using  results  obtained  from  the  solution  of  the  2-D  problem.  The 
present  work  extends  a  simple  yet  accurate  model  developed  recently  for  composite  plates 
and  shells,  namely  Variational  Asymptotic  Plate  and  Shell  Analysis  (VAPAS)  by  Yu  et  al 
(2002a, b,  2003)  so  that  thermoelastic  effects  can  be  treated  in  the  same  framework. 

Since  the  procedure  is  quite  similar,  the  authors  have  chosen  to  repeat  some  formulae 
and  text  from  their  previous  publications  in  order  to  make  the  present  paper  more  self- 
contiiined.  The  present  theory  has  been  implemented  into  the  computer  program  VAPAS. 
The  hygro  effects  due  to  moisture  can  be  treated  in  a  similar  manner  as  thermal  effects. 
Thus,  for  simplicity  of  presentation,  the  hygro  effect  is  not  included  in  the  formulation.  It 
has.  however,  also  been  implemented  in  VAPAS.  Now,  one  can  use  VAPAS  along  some  2-D 
plate  solver  (say,  some  finite  element  program  such  as  DYMORE,  Bauchau  (1998))  to  carry 
out  an  accurate  and  efficient  hygrothermo elastic  analysis  for  composite  plates. 

3-D  Formulation 

A  point  in  the  plate  can  be  described  by  its  Cartesian  coordinates  Xi  (see  Fig.  1),  where 
Xa  are  two  orthogonal  lines  in  the  reference  plane  and  X3  is  the  normal  coordinate.  (Here  and 
throughout  the  paper,  Greek  indices  assume  values  1  and  2  while  Latin  indices  assume  1, 
2,  and  3.  Repeated  indices  are  summed  over  their  range  except  where  explicitly  indicated.) 
Letting  denote  the  unit  vector  along  Xi  for  the  undeformed  plate,  one  can  then  describe 
the  position  of  any  material  point  in  the  undeformed  configuration  by  its  position  vector  r 
from  a  fixed  point  G,  such  that 

r(xi,X2,X3)  =  r(xi,X2)  d-xsba  (1) 

where  r  is  the  position  vector  from  O  to  the  point  located  by  x^  on  the  reference  surface. 
When  the  reference  surface  of  the  undeformed  plate  coincides  with  its  middle  surface,  it 
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naturally  follows  that 


(f(a;i,a;2,  xs))  =  r(xi,a:2)  (2) 

where  the  angle-brackets  denote  the  definite  integral  through  the  thickness  of  the  plate  and 
will  be  used  throughout  the  paper. 

When  the  plate  deforms,  the  particle  that  had  position  vector  f  in  the  undeformed  state 
now  has  position  vector  R  in  the  deformed  plate.  The  latter  can  be  uniquely  determined  by 
the  deformation  of  the  3-D  body.  Similarly,  another  triad  Bj  is  introduced  for  the  deformed 
configuration.  The  relation  between  Bj  and  b;  can  be  specified  by  an  arbitrarily  large 
rotiition  specified  in  terms  of  the  matrix  of  direction  cosines  C{xi,X2)  so  that 

Cij  =  Bihj  (3) 

subject  to  the  requirement  that  B;  is  coincident  with  b,  when  the  structure  is  undeformed. 
Xow  the  position  vector  R  can  be  represented  as 

R(xi.  xo,  xa)  =  R(xi,  xo)  +  X3B3(xi,  Xo)  +  idix^  xo.  X3)Bi(xi,  xo)  (4) 

where  iv-,  is  the  warping  of  the  normal-line  element.  In  the  present  work,  the  form  of  the 
warping  Wi  is  not  assumed,  as  in  most  plate  theories.  Rather,  these  quantities  are  treated  as 
unknown  3-D  functions  and  will  be  solved  for  later.  Eq.  (4)  is  six  times  redundant  because 
of  the  way  warping  introduced.  Six  constraints  are  needed  to  make  the  formulation  unique. 
The  redundancy  can  be  removed  by  choosing  appropriate  definitions  of  R  and  Bj.  One  can 
define  R  similarly  as  Eq.  (2)  to  be  the  average  position  through  the  thickness,  from  which 
it  follows  that  the  warping  functions  must  satisfy  the  three  constraints 


(«;i(xi,X2,X3))  =  0  (5) 

Another  two  constraints  can  be  specified  by  taking  B3  as  the  normal  to  the  reference  surface 
of  the  deformed  plate.  It  should  be  noted  that  this  choice  has  nothing  to  do  with  the 
famous  Kirchhoff  hypothesis.  Indeed,  it  is  only  for  convenience  in  the  derivation.  In  the 
Kirchhoff  assumption,  no  local  deformation  of  the  transverse  normal  is  allowed.  On  the 
other  hand,  according  to  the  present  scheme  we  allow  all  possible  deformation,  classifying 
all  deformation  other  than  that  of  classical  lamination  theory  (CLT)  as  warping,  which  is 
assumed  to  be  small.  This  assumption  is  valid  if  the  strain  is  small  and  if  the  order  of  the 
local  rotation  (i.e.  the  rotation  of  the  normal  line  due  to  warping)  is  of  the  order  of  the  strain 
or  smaller;  see  Danielson  (1991). 

Based  on  the  concept  of  decomposition  of  rotation  tensor,  Danielson  and  Hodges  (1987) 
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and  Danielson  (1991),  the  Jauman-Biot-Cauchy  strain  components  for  small  local  rotation 
are  given  by  , 

Ty  =  \{Fii  +  Fji)  -  S,j  (6) 

where  Fy  is  the  mixed-basis  component  of  the  deformation  gradient  tensor  such  that 

Fij  =  Bi  •  G,g*^  •  b,-  (7) 

Here  G*,  =  ^  is  the  covariant  basis  vector  of  the  deformed  configuration  and  the 
contravariant  base  vector  of  the  undeformed  configuration  and  g^  —  gk  =  bfc.  One  can 
obtain  G^.  with  the  help  of  the  definition  of  so-called  generalized  2-D  strains  similarly  as 
Hodges  et  al.  (1993),  given  by 

IHq;  I 

=  {-Ka3B,3  X  Bs  +  IQzB-i)  x  B,-  (8) 

where  Ea3  and  are  the  2-D  generalized  strains  and  (  ).q  =  Here  one  is  free  to  set 
-12  =  -21? 

Bi  .•  R, 2  =  B2  •  R 1  (9) 

which  can  serve  as  another  constraint  to  specify  the  deformed  configuration. 

With  the  assumption  that  the  strain  is  small  compared  to  unity,  which  has  the  effect  of 
removing  all  the  terms  that  are  products  of  the  warping  and  the  generalized  strains,  one  can 
express  the  3-D  strain  field  as 

r  =  -f-  r,e  -t-  r,,w,i  +  ri,w,2  (lo) 

where 

r  =  [rii  2ri2  r22  2ri3  2r23  rsaj""  (ii) 

W  =  [wi  W2  (12) 

e  =  [sii  2^12  £^22  Fii  K12  +  K21  ^22]^  (13) 
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and  all  the  operators  are  defined  as: 


r. 


0 

0 

0 

a 

dxz 

0 

0 


0 

0 

0 

0 

a 

dxz 

0 


0 

0 

0 

0 

0 

a 

dxz. 


1  0  0  X3  0  0 

0  1  0  0  2:3  0 

0  0  1  0  0  X3 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 


r/x  = 


0  0 
1  0 
0  0 
0  1 
0  0 
0  0 

0  o' 
0  0 
1  0 
0  0 
0  1 
0  0 


(14) 


In  the  present  work,  we  only  consider  one-way  therniomechanical  coupling  and  tempera¬ 
ture  change  due  to  the  deformation  of  the  plate  is  negligible.  Then  we  can  use  the  Helmholtz 
free-energy  functional  (Reddy  (1997))  without  quadratic  terms  involving  temperature  only 
to  carry  out  the  analj-sis.  The  free  energy  per  unit  area  (which  is  the  same  as  the  free  energy 
of  the  normal-line  element)  can  be  written  as 

^  D  r  -  r^Da  (is) 

where  T  is  the  difference  of  temperature  inside  the  structure  with  respect  to  the  reference 
temperature  when  the  plate  is  stress  free,  and  D  is  the  3-D  6x6  material  matrix,  which 
consists  of  elements  of  the  elasticity  tensor  expressed  in  the  global  coordinate  system  Xi.  a 
is  a  6  X  1  column  matrix  representing  the  3-D  thermal  expansion  coefficients.  These  matrices 
are  in  general  fully  populated.  However,  if  it  is  desired  to  model  laminated  composite  plates 
in  which  each  lamina  exhibits  a  monoclinic  sjunmetry  about  its  own  mid- plane  and  is  rotated 
about  the  local  norrrial  to  be  a  layer  in  the  composite  laminated  plate,  then  as  shown  in  Yu 
et  al.  (2002a),  some  parts  of  these  matrices  will  always  vanish  no  matter  what  the  layup 
angle  is. 

To  deal  with  applied  mechanical  loads  as  well  with  thermal  loads,  we  will  at  first  leave 
open  the  existence  of  a  potential  energy  and  develop  instead  the  virtual  work  of  the  applied 
mechanical  loads.  The  virtual  displacement  is  taken  as  the  Lagrangeau  variation  of  the 
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displacement  field,  such  that 


d'R  =  SqgiBi  +  x^SibgiBi  x  B3  +  d'lOjBi  +  x  WjBj  (16) 

where  the  virtual  displacement  of  the  reference  surface  is  given  by 

5?B,=iuB,  (17) 

and  the  virtual  rotation  of  the  reference  surface  is  defined  such  that 

5Bj  =  ^  (18) 

Since  the  strain  is  small,  one  may  safely  ignore  products  of  the  warping  and  the  loading  in 
the  virtual  rotation  term.  Then,  the  work  done  through  a  virtual  displacement  due  to  the 
applied  loads  TiB,-  at  the  top  surface  and  ,3iBi  at  the  bottom  surface  and  body  force  piBi 
through  the  thickness  is 

h 

d'lr  =  {Ti  +  3i  +  {6i))SqBi  +  5i’Ba  ^  ("^o  -  Pa)  +  (XsCpa)  +  S  {TiWf'  +  3iW~  +  (piU?i))  (19) 

where  r,,  Pi,  and  p;  are  taken  to  be  independent  of  the  deformation,  and 

( )”  =  ( )|j,3=_h .  By  introducing  column  matrices  5q,  Sip,  r,  p,  and  d>,  which  are  formed  by 
stacking  the  three  elements  associated  with  indexed  symbols  of  the  same  names,  and  using 
Ecis.  (1),  (3),  and  (4),  one  may  write  the  virtual  work  in  matrix  form,  so  that 

6W  =  5(^ f  +  +  5  {t^w^  +  p^w~  +  (20) 

where 

/  =  r  +  /3+(0) 

'  |(ri  - /3i)  +  (a;30i)  1  ^21) 

m=<  ^{t2  -  P2)  3- {xz(t>2)  > 

0 

The  complete  statement  of  the  problem  can  now  be  presented  in  terms  of  the  principle 
of  virtual  work,  such  that 

5U  -W  =  0  (22) 

In  spite  of  the  possibility  of  accounting  for  nonconservative  forces  in  the  above,  the  problem 
that  governs  the  warping  is  conservative.  Thus,  one  can  pose  the  problem  that  governs  the 
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warping  as  the  minimization  of  a  total  potential  functional 

X[  =  U  +  W  (23) 

so  that 

511  =  0  (24) 

in  which  only  the  warping  displacement  is  varied,  subject  to  the  constraints  Eq.  (5).  This  im¬ 
plies  that  the  potential  of  the  applied  mechanical  loads  for  the  functional  governing  warping 
is  given  by 

W  =  —  0^w~  —  {(jP'w)  (25) 

Below,  for  simplicity  of  terminology,  we  will  refer  to  11  as  the  total  potential  energy,  or  the 
total  energy. 

By  principle  of  minimum  total  potential  energy,  one  can  solve  for  the  unknown  warping 
functions  by  minimizing  the  functional  in  Eq.  (23)  subject  to  the  constraints  of  Ecp  (5). 
Up  to  this  point,  this  is  simply  an  alternative  formulation  of  the  original  3-D  elasticity 
problem.  If  we  attempt  to  solve  this  problem  directly,  we  will  meet  the  same  difficulty  as 
solving  any  full  3-D  elasticity  problem.  Fortunately,  as  shown  in  Yu  et  al.  (2002a, b,  2003), 
the  VAM  can  be  used  to  calculate  the  3-D  warping  functions  asymptotically.  The  through- 
the-thickness  analysis  is  one  dimensional  and  can  be  solved  analjiiically.  However,  finite 
element  discretization  is  preferred  to  solve  the  minimization  problem  for  the  sake  of  dealing 
with  multiple  layers  and  arbitrary  monoclinic  material.  A  5-noded  isoparametric  element  is 
used  because  we  need  the  second-order  warping  functions,  which  are  piecewise,  fourth-degree 
polynomials.  Discretizing  the  transverse  normal  line  into  1-D  finite  elements,  one  can  express 
the  warping  field  as 

w{xi)  =  S{xz)V{xuX2)  (26) 

where  5  is  the  shape  function  and  V  is  the  nodal  value  of  warping  field  along  the  transverse 
normal.  Substituting  Eq.  (26)  into  Eq.  (23),  one  can  express  the  total  energy  in  discretized 
form  as 

2n=V'^EV +  2V'^{Duee  +  DkiXi  +  Dhi^y2) 

+  2(V;[Di,,e-t-VjDi,,e  +  VjDi,i,V2) 
-2VW-2e^ae-2V'[ai,-2Vjai,+2V'^L 
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where  L  contains  the  load  related  terms  such  that 


L  =  -  S-^p  - 


(28) 


The  new  matrix  variables  carry  the  properties  of  both  the  geometry  and  material: 


E  =  ([r,5]^D[rfc5]) 
Dhi,  =  <[r;.5]^D[r,,5]> 
A,  =  {rjDT,) 

Ai.  =  {[ri.sfDTe) 
ah^iiThSfDaT) 
ai,  =  {[ri,SfDaT) 


Dh,  =  {[rHSfDT,) 

Dhi,  =  ([r,5]^D[r,,5]) 

=  {[Ti.Sf  D[Ti,S]) 

Di,i,  =  {nsf  D[ri,S]) 

Dl,e  =  <[n,5]^T>r,> 

c.,  =  (r[Dar> 

oci,  =  {[Ti,SYDaT)  (29) 


Although  the  theory  itself  allows  for  an  analytical  representation  for  arbitrary  temperature 
distribution  through  the  thickness,  here  a  fourth-degree  polynomial  is  used  to  approximate 
the  temperature  distribution  for  each  normal-line  element.  The  discretized  form  of  Eq.  (5) 
is 

V'^Hip  =  0  (30) 

where  H  =  (5^  S)  and  t/’  is  the  normalized  kernel  matrix  of  E  such  that  =  I.  Now 

our  problem  is  transformed  to  minimize  Eq.  (27)  numerically,  subject  to  the  constraints  in 
Eq.  (30). 


Dimensional  Reduction 

To  rigorously  reduce  the  original  3-D  problem  to  a  2-D  plate  problem,  one  must  attempt 
to  reproduce  the  energy  stored  in  the  3-D  structure  in  a  2-D  formulation.  This  dimensional 
reduction  can  only  be  done  approximately,  and  one  way  to  do  it  is  by  taking  advantage 
of  the  smallness  of  h/l.  The  small  pararneter  e,  representing  the  order  of  the  generalized 
2-D  strains  e  has  already  been  taken  advantage  of  when  we  derive  Eq.  (10).  To  reduce  the 
number  of  small  parameters  in  the  asymptotic  analysis,  it  is  reasonable  to  assume  that  the 
order  of  strains  due  to  thermal  loading  is  of  the  order  s.  Thus,  the  quantities  of  interest 
assume  the  following  orders: 

€ai3  ~  flKaP  E  uT  ^  E  fz  ^J■{h/lf  E  ~  f.c(h/l)E  TUa  ~  f.lh{h/l)E  (31) 
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where  i-i  is  the  order  of  the  material  constants  (all  of  which  are  assumed  to  be  of  the  same 
order).  It  is  noted  that  m3  =  0. 

Having  assessed  the  orders  of  all  the  interested  quantities,  the  VAM  can  be  used  to 
mathematically  perform  a  dimensional  reduction  of  the  3-D  problem  to  a  series  of  2-D  models, 
similarly  as  what  have  been  done  in  Yu  et  al  (2002a, b,  2003). 

The  VAM  requires  one  to  find  the  leading  terms  of  the  functional  according  to  the  different 
orders.  Since  only  the  warping  is  varied,  the  leading  terms  needed  are  all  of  those  terms 
associated  with  warping.  For  the  zeroth-order  approximation,  these  leading  terms  of  Eq.  (27) 
are 

2U*o  =  V^EV +2V^Dhee-2V'^ah  (32) 

The  Euler-Lagrange  equation  for  functional  Eq.  (32)  subject  to  constraints  Eq.  (30)  can  be 
obtained  by  usual  procedure  of  calculus  of  variation  with  the  aid  of  a  Lagrange  multiplier  as 
follows: 

EV  +  Dhee-ah  =  Hi^u\  (33) 

Considering  the  properties  of  the  kernel  matrix  0,  one  calculates  the  Lagrange  multiplier  A 
as 

A  =  ip'^{Dhe€  -  ah)  (34) 

Substituting  Eq.  (34)  back  iirto  Eq.  (33),  we  obtain 

EV  =  {Hv)V  -I){Dhee-ah)  (35) 

There  is  a  unique  solution  of  zeroth-order  warping  functions  and  can  be  written  as: 

V  =  VQe  +  VT  =  Vo  (36) 

Substituting  Eq.  (36)  back  into  Eq.  (27),  one  can  obtain  the  total  energy  asymptotically 
correct  up  to  the  order  of  as 

2R^  =  t^(V^D^  +  D„)e  +  t'^(Dl,VT-VoaH-‘ia,)  (37) 

Xote  the  quadratic  terms  associated  with  temperature  -V^ah  is  dropped  due  to  the  same 
reason  for  Eq.  (15).  This  2-D  free  energy,  Eq.  (37),  is  the  same  as  what  is  used  in  CLT  for 
thermoelastic  analysis: 

2110  =  -  2€^iV-r  (38) 
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with 


.4  =  {Vq  Dhe  +  Ae)  =  «£  +  -^{Vo  Oih  -  DI^Vt)  (39) 

Although  the  energy  of  this  approximation  coincides  with  CUT,  we  have  not  used  any  ad  hoc 
kinematic  assumptions  such  as  the  KirchhofF  assumption  to  obtain  this  result.  Moreover, 
the  transverse  normal  strain  from  our  zeroth-order  approximation  is  not  zero. 

It  is  understood  that  our  zeroth-order  approximation  will  give  the  same  Stress  results 
as  what  is  obtained  from  CLT,  i.e.,  all  the  transverse  stress  components  which  are  very 
important  for  analyzing  the  failure  of  composite  plates  cannot  be  predicted.  One  must  carry 
out  the  next  approximation  so  that  those  quantities  can  be  approximately  predicted.  To 
obtain  the  first-order  approximation,  we  simply  perturb  the  zeroth-order  result,  resulting  in 
warping  functions  of  the  form 

y  =  V^o  +  Ki  (40) 

Substituting  Eq.  (40)  back  into  Eq.  (10)  and  then  into  Eq.  (27),  one  can  obtain  the  leading 
terms  for  the  first-order  approximation  as 

2ni  =  EVi  +2V^  Die, I  +  2V[  Ae  2  +  2V^Lt  +  2  L  (41) 

with 

Di  =  {Dm,-DIJVo-Di,, 

Do  =  {Dhi^  —  D]^i^)Vq  —  Aie 

Lt  =  {Dhh  -  DliJVr,i  +  {Dhi2  -  DIi^)Vt,2  +  +  0^2.2  (42) 

Integration  by  parts  with  respect  to  the  in-plane  coordinates  is  used  here  and  hereafter 
whenever  it  is  convenient  for  the  derivation,  because  the  present  goal  is  to  obtain  an  interior 
solution  for  the  plate  without  consideration  of  boundary  layer  effects. 

Similarly  as  in  the  zeroth-order  approximation,  one  can  solve  the  first-order  warping  field 
as 

Vi  =  +  Vl2^,2  +  +  ^IL  (43) 

and  obtain  a  total  energy  that  is  asymptotically  correct  up  to  the  order  of  i.i{h/iy^£,  given 
by 

2ni  =  e'^Ae  4-  1  +  2  +  eJPe.o  +  2e^F  -  2e^FT  (44) 
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where 


B  =  !/„’' All, Vi + 


V  =  VjDi,,,Vo  +  VlDi 


C  =  Vi  D,,,,Vo  +  i(ViD,  +  D\  Ki.) 

F  =ViL  -  i(OrVu,,  +  V^i,,  +  +  V,|t,2)  (45) 


with  the  non-mechanical  load  due  to  temperature 


Ft  —Nt  +  Vq  AwiVr.n  +  Vj A2i2^,22  +  Vq  +  Df^i^)VT,i2 

+  2(^n^T,l  +  y\2^T2  +  D^ViT,l  +  D2  ViT,2) 


(46) 


Here  the  monoclinic  symmetry  has  already  been  taken  advantage  of  to  obtain  the  asymptot¬ 
ically  correct  energy  in  Eq.  (44).  The  applied  mechanical  loads  and  temperature  distribution 
should  not  vary  rapidly  over  the  plate  surface;  otherwise  the  structure,  although  plate-like, 
can  not  be  analyzed  with  enough  accuracy  using  a  reduced  plate  model. 


Transforming  Into  Reissner-Mindlin  Model 
Although  Eq.  (44)  is  asymptotically  correct  through  the  second  order  and  straightfor¬ 
ward  use  of  this  free  energy  expression  is  possible,  it  involves  more  complicated  boundary 
conditions  than  necessary  since  it  contains  derivatives  of  the  generalized  strain  measures. 
To  obtain  an  energy  functional  that  is  of  practical  use,  one  can  transform  the  present  ap¬ 
proximation  into  a  Reissner-Mindlin  type  model.  We  would  like  to  state  that  fitting  the 
asymptotic  energy  into  such  model  is  just  a  choice,  and  the  possibility  of  fitting  the  same 
energy  into  other  more  sophisticated  plate  models  is  under  investigation. 

In  a  Reissner-Mindlin  model,  there  are  two  additional  degrees  of  freedom,  which  are  the 
transverse  shear  strains.  These  are  incorporated  into  the  rotation  of  transverse  normal.  If 
we  introduce  another  triad  B*  for  the  deformed  Reissner-Mindlin  plate,  the  definition  of  2-D 
strains  becomes 

R,  =  B:-h£::^B^  +  27a3B^ 

X  BJ  +  X  B*  (47) 

where  the  transverse  shear  strains  are  7  =  [2713  2723]^-  From  the  definitions  in  Eqs.  (8) 
and  (47),  one  can  obtain  the  Rodrigues  parameters  corresponding  to  the  rotation  relating 
Bi  and  B^.  Using  the  procedures  listed  in  Hodges  (1987),  one  can  express  the  classical 
strain  measures  e  in  terms  of  the  strain  measures  of  the  Reissner-Mindlin  plate  model  (see 
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Yu  (2002)  for  the  details  of  the  derivation) : 

€  —  TZ  —  T)ci''l\a 


(48) 


where 


tT 


= 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

_ 

T 

^2  = 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

n  =  u-:, 

2, 

-12 

-22 

Kl,+Kl, 

(49) 


Now  one  can  express  the  energy,  Eq.  (44),  correct  to  second  order,  in  terms  of  strains  of  the 
Reissner-Mindlin  model  as 


2ni  =71^ ATI  -  +  2nlcn2  +  7^5^7^,2 

+  2n^F  -  2R^Ft  +  2-iiDaNT 


(50) 


The  generalized  Reissner-Mindlin  model  used  in  the  2-D  thermoelastic  analyses  is  of  the 


form 


2^7^  ^  FJ'A'R  7^C?7  +  271^ Fti  +  2-i^Fy  +  271^ Fjn  +  2-/^ Ft-, 


(51) 


To  find  an  equivalent  Reissner  model  Eq.  (51)  for  Eq.  (50),  one  has  to  eliminate  all  partial 
derivatives  of  the  classical  2-D  strain  measures.  The  equilibrium  equations  are  used  to 
achieve  this  purpose.  Prom  the  two  equilibrium  equations  balancing  bending  moments  with 
applied  moments  which  is  calculated  from  Eq.  (21),  one  can  obtain  the  following  formula 


G'y  +  F-f  +  Ft-,  =  7?^{A7Z^a  +  FTi,a  +  Fth^oc)  + 
Using  Eq.  (52),  one  can  rewrite  Eq.  (50)  as 


{:} 


(52) 


where 


2n,  =  n^ATl  +  ■fG't  +  2R^F  -  2VJFt  -  2-/V^Nr.^  +  U’ 


U‘ =  11^1311,^+211^10113  +  1^2^11.2 


(53) 


(54) 
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and 


B  =  B  +  AViG-^VfA 

C  =  C  +  AViG-^VlA 

D  =  D  +  AV2G-^VIA  (55) 

If  we  can  drive  U*  to  zero  for  any  TZ,  then  we  have  found  an  asymptotically  correct  Reissner- 
Mindlin  plate  rhodel.  For  generally  anisotropic  plates,  this  term  will  not  be  zero;  but  we 
can  minimize  the  error  to  obtain  a  Reissner-Mindlin  model  that  is  as  close  to  asymptotical 
correctness  as  possible.  The  accuracy  of  the  Reissner-Mindlin  model  depends  on  how  close 
to  zero  one  can  drive  this  term  of  the  energy. 

One  could  proceed  with  the  optimization  at  this  point,  but  the  problem  will  require  a 
least  scpiares  solution  for  3  unknowns  (the  shear  stiffness  matrix  G)  from  a  linear  system  of  78 
eciuations  (12x12  and  symmetric).  This  optimization  problem  is  too  rigid.  The  solution  will 
be  better  if  we  can  bring  more  unknowns  into  the  problem.  As  stated  in  Sutyrin  and  Hodges 
(1996),  there  is  no  unique  plate  theory  of  a  given  order.  One  can  relax  the  constraints  in 
Eq.  (5)  to  be  (ic,)  =  const  and  still  obtain  an  asymptotically  correct  strain  energy.  Since  the 
zeroth-order  approximation  gives  us  an  asymptotic  model  corresponding  to  classical  plate 
theory,  we  only  relax  the  constraints  for  the  first-order  approximation.  This  relaxation  will 
modify  the  warping  field  to  be 

1/ 1  =  -l-  Vj[2C,2  ViL  +  I'lr  +  Lic^i  -1-  ^2^,2  (56) 

where  Ii,  L2  consist  of  24  constants.  The  remaining  energy  U*  will  also  be  modified  to  be 

C/*  =  -f  27^^C7^,2  +  7^5i)7^,2  (57) 

and 

B  =  B  +  2LjDi  C  =  C  +  {L(D2  +  Di'^L2)  D  =  D  +  2^D2  (58) 

Since  now  we  have  27  unknowns,  the  optimization  is  much  more  flexible.  It  can  give  us  a 
more  optimal  solution  for  the  shear  stiffness  matrix  G  to  fit  the  second-order,  asymptoti¬ 
cally  correct  energy  into  a  Reissner-Mindlin  model.  In  other  words,  here  we  have  found  the 
Reissner-Mindlin  model  that  describes  as  closely  as  possible  the  2-D  energy  that  is  asymptot¬ 
ically  correct  through  the  second  order  in  h/l.  However,  the  asymptotical  correctness  of  the 
warping  field  to  that  same  order  can  only  be  ascertained  after  obtaining  another  higher-order 
approximation,  which  will  be  discussed  in  the  next  section. 
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And  after  minimizing  U*.  the  "best”  total  energy  to  be  used  for  the  2-D  plate  Reissner- 
Mindlin  model  can  be  expressed  as: 

2^7^  =  7^^A7^  +  7^G7  +  27^^F-27^^F^-27^Pa^T,a  (59) 

Recovering  Relations 

From  the  above,  we  have  obtained  a  Reissner-Mindlin  plate  model  which  is  as  close  as  pos¬ 
sible  to  being  asymptotically  correct  in  the  sense  of  matching  the  total  energy.  The  stiffness 
matrices  A,  G,  load-related  terms,  and  non- mechanical  stress  resultants  can  be  used  as  input 
for  a  plate  theory  derived  from  the  total  energy  obtained  here.  The  geometrically  nonlinear 
theory  presented  in  Hodges  et  al.  (1993)  is  an  appropriate  choice,  but  some  straightforward 
modification  of  the  loading  terms  is  required. 

However,  while  it  is  necessary  to  accurately  calculate  the  2-D  displacement  field  of 
composite  plates,  this  is  not  sufficient  in  many  applications.  Ultimately,  the  fidelity  of 
a  reduced-order  model  such  as  this  depends  on  how  well  it  can  predict  the  3-D  results  in 
the  original  3-D  problem.  Hence,  recovering  relations  should  be  provided  to  complete  the 
reduced-order  model.  By  recovering  relations,  we  mean  expressions  for  3-D  displacement, 
strain,  and  stress  fields  in  terms  of  2-D  quairtities  and  0:3.  For  validation,  results  obtained 
for  the  3-D  field  variables  from  the  reduced-order  model  must  be  compared  with  those  of  the 
original  3-D  model. 

For  an  energy  that  is  asymptotically  correct  through  the  second  order,  we  can  recover  the 
3-D  displacement,  strain  and  stress  fields  only  through  the  first  order  in  the  strict  sense  of 
asymptotical  correctness.  Using  Eqs.  (1),  (3),  and  (4),  one  can  recover  the  3-D  displacement 
field  through  the  first  order  as 

G31 

Uzd  =  tt2d  -h  X3  Czz  +  SVq  -b  SV 1  (60) 

Czz  -‘1. 

where  Uzd  is  the  column  matrix  of  3-D  displacements  and  U2d  is  the  plate  displacements. 
Cij  are  the  components  of  global  rotation  tensor  from  Eq.  (3).  And  from  Eq.  (10),  one  can 
recover  the  3-D  strain  field  through  the  first  order  as 

r  =  r,,5(Uo  +  Ui)  +  r,e  +  r,^5Uo,i  +  r,35Uo,2  (6i) 
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Then,  one  can  use  the  inverse  of  3-D  Duhamel-Neuniann  law 

a  =  DT-DaT  (62) 

to  obtain  3-D  stresses  aij. 

Since  we  have  obtained  an  optimurn  shear  stiffness  matrix  G,  some  of  the  recovered  3-D 
results  through  the  first  order  are  better  than  classical  theory  and  conventional  first-order 
deformation  theory.  However,  for  the  transverse  normal  component  of  strain  and  stress  {i.e. 
r33  and  (T33),  the  agreement  is  not  satisfactory  at  all.  Let  us  recall,  that  the  Reissner-Mindlin 
theory  that  has  been  constructed  only  ensures  a  good  fit  with  the  asymptotically  correct  3- 
D  displacement  field  of  the  first  order  (while  energy  is  approximated  to  the  second  order). 
Thus,  in  order  to  obtain  recovering  relations  that  are  valid  to  the  same  order  as  the  energj'', 
the  V'4M  iteration  needs  to  be  applied  one  more  time. 

Using  the  same  procedure  listed  in  previous  section,  the  second-order  warping  can  be 
obtained  and  expressed  symbolically  as 

V2  =  Voie.li  +  V22C.12  +  ^^3^.22  +  ^2L  +  (^3) 

Eq.  (63)  is  obtained  by  taking  the  original  first-order  warping  Vi  to  be  the  result  of  the 
first-order  approximation.  It  is  clear  that  V2  is  one  order  higher  than  Vi  which  confirms 
that  Vi  is  the  first-order  approximation.  One  might  be  tempted  to  use  Vi  in  the  recovering 
relations.  However,  the  VAM  has  split  the  original  3-D  problem  into  two  sets  of  problems. 
As  far  as  recovering  relations  concerned,  it  is  observed  that  the  normal-line  analysis  can 
at  best  give  us  an  approximate  shape  of  the  distribution  of  3-D  results.  The  2-D  plate 
analysis  will  govern  the  global  behavior  of  the  structure.  Since  V 1  is  the  warping  that  yields 
a  Reissner-Mindlin  plate  model  that  is  as  close  as  possible  to  being  asymptotically  correct, 
we  should  still  use  Fi  in  the  recovering  relations  instead  of  Vi.  By  doing  this,  we  choose  to 
be  more  consisterit  with  Reissner-Mindlin  plate  model  while  compromising  somewhat  on  the 
asymptotical  correctness  of  the  shape  of  the  distribution.  It  has  been  verified  by  numerical 
examples  that  this  is  a  good  choice. 

Hence,  we  write  the  3-D  recovering  relations  for  displacement  through  the  second  order 
as  ■ 

'  G31  ]  _ 

Uid  —  U2d  +  ^3  *  G32  /  +  5'(Vo  +  Vi  +  V2)  (64) 

,G33-iJ 
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and  the  strain  field  through  the  second  order  is 

r  =  ThS{Vo  +  Vi+  V2)  +  r^e  +  r[,S{VoA  +  Vi,i)  +  Ti,S{Vo,2  +  ^1,2)  (65) 

Again  the  stresses  through  the  second  order  can  be  obtained  from  use  of  the  3-D  constitutive 
law,  Eq.  (62). 

Numerical  Examples 

The  computer  program  VAPAS  has  been  extended  to  implement  the  present  theory. 
Several  numerical  examples  are  given  here  to  validate  the  proposed  theory  and  code  against 
the  3-D  thermoelasticity  solutions. 

First  to  assess  the  asymptotical  correctness  of  the  proposed  theory,  we  study  a  cylindrical 
bending  type  problem  for  an  isotropic  plate  with  E  as  the  Young’s  modulus,  i/  Poisson’s 
ratio  and  a  thermal  expansion  coefficient.  The  plate  is  simply  supported  with  width  L  along 
;ri  axis  (the  "lateral”  direction)  and  infinitely  long  along  the  X2  axis  (the  ‘longitudinal” 
direction)  under  the  following  temperature  changes: 

T  =  Toixz)  sin  (66) 

Let  us  assume  the  thickness  of  the  plate  is  h,  and  the  normalized  thickness  coordinate 
C  =  Xi/h,  then  the  small  parameter  used  in  our  theory  is: 

'’  =  7  =  T 

When  there  is  a  uniform  temperature  Tc  change  through  the  thickness,  the  nontrivial  dis¬ 
placements  and  stresses  are  listed  in  Table  1.  The  exact  solutions  are  obtained  based  on 
Savoia  and  Reddy  (1995)  and  expanded  into  a  series  in  terms  of  5  with  o(*)  denoting  terms 
asymptotically  smaller  than  the  order  of  *.  The  present  theory  can  predict  the  correct  re¬ 
sults  up  to  the  second  order  of  S  with  respect  to  the  leading  terms  for  each  3-D  quantities, 
which  clearly  demonstrate  that  our  theory  is  asymptotically  correct  up  to  the  second  order 
for  this  particular  problem.  We  admit  that  the  prediction  of  transverse  components  for  this 
problem  is  out  of  the  power  of  our  theory.  However,  this  should  not  mislead  the  reader  to 
assume  that  the  present  theory  is  asymptotically  correct  up  to  the  second  order  for  general 
cases.  The  authors  are  aware  that  the  proposed  theory  can  be  at  best  asymptotically  correct 
up  to  the  second  order  for  particular  cases.  For  general  cases,  however,  the  theory  can  only 
be  interpreted  as  a  Reissner-Mindlin  type  theory  which  is  closest  to  being  asymptotically 
correct.  To  illustrate  the  above  statement,  we  provide  the  results  for  the  same  isotropic  plate 
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under  a  linearly  distributed  temperature  through  the  thickness,  assuming  To  =  QTi.  Here 
only  nontrivial  displacement  results  are  listed  in  Table  2  which  is  sufficient  for  the  aforemen¬ 
tioned  purpose.  One  can  observe  from  Table  2  that  there  is  a  slight  difference  for  the  second 
order  prediction  between  the  present  theory  and  exact  solution.  It  is  interesting  to  note  that 
if  one  sets  i/  equal  to  zero  the  difference  disappears.  Evidently  some  information  belonging 
to  second  order  and  indeed  included  in  the  asymptotically  correct  energy  cannot  be  captured 
in  a  Reissner-Mindlin  type  model.  When  we  transform  the  asymptotically  correct  model  into 
a  Reissner-Mindlin  model,  this  information  is  lost. 

The  present  theory  is  formulated  with  sufficient  generality  to  carry  out  a  thermoelas¬ 
tic  analysis  for  arbitrary  composite  laminated  plates  made  with  monoclinic  material  with 
a  computational  cost  equivalent  to  that  of  first-order  shear-deformation  theory  (FOSDT). 
Hence,  before  overemphasizing  the  loss  of  information  in  the  repackaging  of  the  model  one 
should  determine  if  this  loss  is  exhibited  in  numerical  results.  To  investigate  this,  we  will 
present  some  numerical  results  for  composite  plates  to  demonstrate  the  advantages  of  our 
theory  relative  to  CLT  and  FOSDT.  The  plate  we  are  going  to  study  has  material  properties 
given  by 


El  =  172.4  GPa  (25  x  10®  psi) 
Glt  =  3.447  GPa  (0.5  x  10®  psi) 
~  0.25 

=  0.139  X  10“®/ “C 


Et  =  6.895  GPa  (10®  psi) 

Gtt  =  1-379  GPa  (0.2  x  10®  psi) 
utt  ”  0.25 

ar  =  9  X  10"®/  (68) 


For  the  purpose  of  comparing  with  the  exact  solution,  we  still  consider  cylindrical  bending 
U'pe  problem.  In  lieu  of  our  definition  of  small  parameter  Eq.  (67),  even  if  our  theory  is 
asjanptotically  correct  up  to  the  second  order,  we  require  6^  «  1.  If  we  assume  the  thickness 
is  25.4  mm  and  L  =  254  mm,  5  will  be  approximately  0.314.  which  can  be  considered  small 
in  our  theory.  Two  different  cases  are  investigated: 

•  case  1:  nearly  cross  ply,  [—0.5°/89.5°]  under  To  =  Tc +  C^c 

•  case  2:  nearly  cross  ply,  [90.5°/0.5°/90.5°/0.5°]  under  To  =  Tc  +  -}- C^Tc -1- C^Tc 

T3  =  03  =  ^  sin  (^) ,  Ta  =  Pa  =  0,  and  po  =  EtoltTc/^. 

Because  thermal  stresses  due  to  temperature  change  are  the  most  interesting  quantities, 
we  only  present  stress  results  here  with  mentioning  that  the  accuracy  of  displacements  and 
strains  is  similar  to  that  of  stresses.  For  ease  1,  results  from  VAPAS  (dots  in  the  plots), 
are  compared  with  those  from  CLT  (dash-dotted  line),  FOSDT  (dashed  line)  and  the  exact 
solution  (solid  line)  in  Figs.  2-7.  Note  that,  because  the  2-D  variables  are  either  sine 
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functions  or  cosine  functions,  cToj  and  are  plotted  at  Xi  =  Ljl  and  cTos  at  Xi  =  L.  The 
results  presented  here  are  normalized  as  follows: 


9aij 

EtCHtTc 


(69) 


As  one  can  observe  from  the  plots,  for  VAPAS  results  are  much  closer  to  the  exact 
solutions  than  CLT  and  FOSDT.  VAPAS  also  does  a  fairly  good  job  for  predicting  the 
transverse  stress  components  which  for  the  isotropic  plate  under  uniform  temperature  change 
we  concluded  was  out  of  the  power  of  VAPAS  because  these  terms  are  asymptotically  smaller 
than  the  second  order.  One  can  infer  that  due  to  the  special  layup  scheme  (cross-ply)  the 
dominant  terms  of  trans\^rse  stress  components  could  now  be  asymptotically  eciual  or  larger 
than  the  second  order.  Considering  the  smallness  of  <733  in  comparison  to  the  in-plane 
components,  we  expect  the  slight  shift  of  VAPAS  results  for  this  quantity  to  be  tolerable  for 
most  engineering  applications. 

Finally,  to  demonstrate  that  VAPAS  can  handle  the  temperature  change  through  the 
thickness  exactly  up  to  a  fourth-degree  polynomial  and  both  mechanical  loads  and  thermal 
load  can  be  treated  simultaneously,  we  present  the  results  for  case  2  in  Figs.  8  -  13.  Except 
for  a  small  shift  for  transverse  normal  stress,  all  the  other  results  from  VAPAS  are  almost  on 
top  of  the  exact  solutions.  Careful  readers  may  even  find  there  is  a  small  discontinuity  for  0-33 
which  should  not  be  the  case  in  reality.  The  reason  is  due  to  that  the  stress  results  obtained  by 
VAPAS  are  calculated  directly  using  3-D  constitutive  law  from  the  approximate  strain  field. 
The  approximation  in  the  strain  field  may  cause  the  discontinuity  for  the  transverse  stress 
components.  It  is  worthy  to  emphasize  that  integration  of  the  3-D  equilibrium  ecpiatioiis 
through  the  thickness  is  not  used  here  to  obtain  results  for  the  transverse  stresses  presented 
herein,  in  contrast  to  what  is  usually  done  in  CLT  and  FOSDT. 

Mathematically,  the  accuracy  of  the  present  theory  should  be  comparable  to  that  of 
a  reduced  layer- wise  plate  theory  with  assumed  in-plane  displacements  as  layer- wise  cubic 
polynomials  of  the  thickness  direction  and  transverse  displacement  as  a  layer-wise  fourth- 
degree  polynomial.  However,  the  present  theory  is  still  an  equivalent  single-layer  theory,  and 
the  computational  requirement  is  much  less  than  that  for  layer-wise  theories. 


Conclusion 

A  Reissner-Mindlin  type  plate  model  capable  of  performing  a  thermoelastic  stress  analysis 
of  laminated  composite  plates  has  been  constructed  by  the  variational-asymptotic  method. 
A  general  2-D  constitutive  law  being  as  close  to  asymptotical  correctness  as  possible  has  been 
obtained  by  solving  the  1-D  through-the-thickness  analysis.  The  original  3-D  results  have 
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been  reproduced  as  accurate  as  possible  from  a  Reissner-Mindlin  type  plate  analysis.  The 
resulting  theory  is  as  simple  and  efficient  as  a  first-order  shear  deformation  theory  (FOSDT) 
while  it  also  has  an  accuracy  comparable  to  higher-order  layerwise  theories. 

The  present  study  distinguishes  from  previous  work  reported  in  the  literature  at  least  in 
the  following  four  aspects: 

1.  The  present  formulation  is  in  an  intrinsic  form  which  is  suitable  for  geometrically 
nonlinear  plate  theory  as  well  as  linear  theory. 

2.  The  dimensional  reduction  from  3-D  to  2-D  is  carried  out  systematically  by  using 
variational  asymptotic  method,  w'hich  is  completely  different  from  models  that  rely  on 
the  introduction  of  ad  hoc  kinematic  assumptions  to  reduce  the  dimensionality. 

3.  To  create  a  smooth  interface  with  well-established  2-D  solvers,  the  degrees  of  freedom  of 
the  present  model  are  chosen  to  be  essentially  the  same  as  those  of  traditional  Reissner- 
Mindlin  type  plate  theory.  However,  the  present  model  is  not  a  FOSDT.  The  present 
theory  differs  from  FOSDT  by  representing  all  the  deformations  that  are  purposely 
eliminated  in  the  development  of  FOSDT.  This  is  accomplished  through  allowing  ail 
possible  deformation  in  the  3-D  warping  functions,  which  are  solved  in  turn  by  the 
variational-asymptotic  method. 

4.  The  present  study  has  treated  both  mechanical  and  tliermal  loading.  The  temperature 
distribution  through  the  thickness  can  be  arbitrary  and  is  approximated  in  VAPAS 
by  a  layerwise  fourth-degree  polynomial.  This  is  more  realistic  and  accurate  than 
most  published  models,  in  which  the  temperature  is  assumed  to  be  distributed  linearly 
through  the  thickness  of  the  whole  plate  (single-layer  theories)  or  a  layer  (for  layerwise 
theories). 

The  hygro  effect  due  to  moisture  to  composite  plates  can  also  be  handled  in  exactly  the 
same  procedure  except  one  has  to  replace  the  thermal  expansion  coefficients  with  hygroscopic 
expansion  coefficients  and  temperature  with  moisture.  The  computer  program  VAPAS  can 
now  be  used  along  with  a  2-D  solver  to  perform  an  efficient  yet  accurate  and  detailed  analysis 
for  hygrothermoelastic  behavior  for  laminated  composite  plates  under  severe  environments. 
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Table  1:  3-D  displacements  and  stresses  under  uniform  temperature  change  through 

the  thickness 


Normalized  lateral  displacement  (j;^) 

—  0  ol  0  t"  5760(i/-l)  ^  J 


Exact 

Present 


+  o(S) 


Normalized  transverse  displacement  {j;^) 

r.-  I  nr  I  ^-C(4C^-l)-^(t^+l)j2  I  ^N18C^(6C^-5)(t/+l)^-t^l4t/+l_51  , 

- 2MI^)  °  +  5760(t/-l)a-‘  +O[0  ) 


Exact 

Present 


N ormalized  lateral  in-plane  stress 


Exact 

Present 


;t^(12C— Drf  :r-*(240C-*-120C^+7)  r4  .  ( ):A\ 

■  24(i/-ir®  2880(i/-l)  ^  ) 


Normalized  longitudinal  in-plane  stress 
Exact  -1  - 

Present  -1  -  +  q(.^") 


Normalized  lateral  transverse  shear  stress 


Exact 

Present  Not  available 


N ormalized  transverse  normal  stress  (•^^) 


Exact 

Present  Not  available 
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Table  2:  3-D  displacements  under  linearly  distributed  temperature  change  through 

the  thickness 


Exact 

Present 


Exact 

Present 


Normalized  lateral  displacement  j-^^) 


IT 


C(t'+1)  f 


r_l  _  n-C(20C^-3)(l/+l)  r  .  TtC^i/{\].u-^+2u^-\-'6u^  +  \Av-Z)  r  ,  /  r\ 

_ f _ 120  ^  30(llt/-*-12i/3+34t/^-r2t/+ll)^  ^  ' 

Normalized  transverse  displacement 

+  i(20C-  -  l)(i^  +  1)  +  o(J'^) 


TT 

i^r-2 


A-2  o.  J./’9nr2  _  1  W,y  4.  1  ^  _  t'C  1 1 +2v^ +8t/-  + 14^/-3)  ,  /  v-0\ 


25  OF  39 


List  of  Figures 

1  Schematic  of  plate  deformation  .  27 

2  Distribution  of  the  3-D  stress  an  vs  the  thickness  coordinate  (case  1)  .  .  .  .  28 

3  Distribution  of  the  3-D  stress  ai2  vs  the  thickness  coordinate  (case  1)  .  .  .  .  29 

4  Distribution  of  the  3-D  stress  a22  vs  the  thickness  coordinate  (case  1)  .  .  .  .  30 

5  Distribution  of  the  3-D  stress  ais  vs  the  thickness  coordinate  (case  1)  ....  31 

6  Distribution  of  the  3-D  stress  <723  vs  the  thickness  coordinate  (case  1)  .  .  .  .  32 

7  Distribution  of  the  3-D  stress  <733  vs  the  thickness  coordinate  (case  1)  .  .  .  .  33 

8  Distribution  of  the  3-D  stress  an  vs  the  thickness  coordinate  (case  2)  ...  .  34 

9  Distribution  of  the  3-D  stress  <712  vs  the  thickness  coordinate  (case  2)  ...  .  35 

10  Distribution  of  the  3-D  stress  <722  vs  the  thickness  coordinate  (case  2)  ...  .  36 

11  Distribution  of  the  3-D  stress  C7i3  vs  the  thickness  coordinate  (case  2)  .  .  .  .  37 

12  Distribution  of  the  3-D  stress  <723  vs  the  thickness  coordinate  (case  2)  .  .  .  .  38 

13  Distribution  of  the  3-D  stress  <733  vs  the  thickness  coordinate  (case  2)  ...  .  39 


26  OF  39 


I 


i 


Figure  2: 
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Distribution  of  the  3-D  stress  cn  vs  the  thickness  coordinate  (case  1) 
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In-Plane  Shear  Stress 


Figure  3:  Distribution  of  the  3-D  stress  cri2  vs  the  thickness  coordinate  (case  1) 
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Longitudinal  In-Plane  Stress 


Figure  4;  Distribution  of  the  3-D  stress  <722  vs  the  thickness  coordinate  (case  1) 
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Lateral  Transverse  Shear  Stress 


Figure  5: 


Distribution  of  the  3-D  stress  cris  vs  the  thickness  coordinate  (case  1) 


31  OF  39 


Longitudinal  Transverse  Shear  Stress 


Figure  6: 
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Distribution  of  the  3-D  stress  <723  vs  the  thickness  coordinate  (case  1) 
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Figure  7:  Distribution  of  the  3-D  stress  0^33  vs  the  thickness  coordinate  (case  1) 
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Lateral  In~Plane  Stress 


Figure  8: 
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Distribution  of  the  3-D  stress  crn  vs  the  thickness  coordinate  (case  2) 


34  OF  39 


In-Plane  Shear  Stress 


Figure  9: 


Distribution  of  the  3-D  stress  ai2  vs  the  thickness  coordinate  (case  2) 
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Longitudinal  In-Plane  Stress 


Figure  10:  Distribution  of  the  3-D  stress  G22  vs  the  thickness  coordinate  (case  2) 
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Distribution  of  the  3-D  stress  c^^Q  vs  the  thickness  coord 


Longitudinal  Transverse  Shea 


Normalized  Thickness  Coordinate 

Distribution  of  the  3-D  stress  cr^z  vs  the  thickness  coordinate  (case  2) 


Transverse  Normal  Stress 


Figure  13:  Distribution  of  the  3-D  stress  <733  vs  the  thickness  coordinate  (case  2) 
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A  Simple  Thermopiezoelastic  Model  for  Composite  Plates 
with  Accurate  Stress  Recovery! 

Wenbin  Yu  I  and  Dewey  H.  Hodges§ 

Georgia  Institute  of  Technology,  Atlanta,  Georgia  30332-0150 

E-mail:  dewey . hodges@ae . gatech . edu 

Abstract  A  Reissner-Mindlin  type  model  for  analyzing  laminated  composite  plates 
including  piezoelectric  layers  under  mechanical,  thermal  and  electric  loads  has  been 
constructed  by  the  variational  asymptotic  method.  The  present  work  formulates  the  original 
nonlinear,  three-dimensional,  one-way  coupled,  therniopiezoelasticity  problem  allowing  for 
arbitrary  deformation  of  the  normal  line  and  using  a  set  of  intrinsic  variables  defined  on 
the  reference  plane.  The  variational  asymptotic  method  is  used  to  rigorously  split  the  three- 
dimensional  problem  into  two  problems:  a  nonlinear,  two-dimensional,  plate  analysis  over  the 
reference  plane  to  obtain  the  global  deformation,  and  a  linear  analysis  through  the  thickness 
to  provide  both  the  two-dimensional  generalized  constitutive  model  and  recovering  relations 
I  to  approximate  the  original  three-dimensional  results.  The  obtained  asymptotically  correct 
second-order  free  energy  is  cast  into  the  form  of  the  comnnonly-used  Reissner-Mindlin  type 
model  to  account  for  transverse  shear  deformation.  The  present  theory  is  implemented  into 
the  computer  program  VAPAS  (Variational  Asymptotic  Plate  and  Shell  Analysis).  Results 
for  several  cases  obtained  from  VAPAS  are  compared  with  the  exact  thermopiezoelasticity 
solutions,  classical  lamination  theory  and  first-order  shear-deformation  theory  for  the  purpose 
of  demonstrating  advantages  and  limitations  of  the  proposed  theory.  The  proposed  theory  can 
achieve  an  accuracy  comparable  to  higher-order  layerwise  theories  at  the  cost  of  a  first-order 
shear  deformation  theory. 


Submitted  to:  Smart  Materials  and  Structures 

1.  Introduction 

Research  on  smart  structures  has  received  enormous  attention  in  recent  years  [I,  2,  3,  4,  5]. 
Smart  structures  are  capable  of  sensing  and  reacting  to  external  disturbances  and  thus 
create  the  possibility  of  building  structures  that  are  self-monitorable  and  self-controllable. 
Such  smart  structures  are  promising  candidates  to  meet  the  demanding  requirements  of 
high-strength,  high-stiffness,  and  light-weight  structures  for  modern  engineering,  especially 
aerospace  applications.  Among  many  possible  candidates  for  actuators  and  sensors, 
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piezoelectric  materials  receive  the  most  attention.  One  reason  for  this  preference  is  that 
piezoelectric  materials  can  directly  relate  electrical  signals  to  strains  in  the  material  and  vice 
versa.  Thus,  they  can  be  used  both  as  actuators  and  sensors.  Moreover,  in  most  cases, 
piezoelectric  materials  are  used  along  with  tailored  anisotropic  materials  to  maximize  the 
intelligence  of  a  smart  structure.  While  most  research  has  been  focused  on  the  behavior  of 
piezoelectric  structures  under  isothermal  conditions,  an  increasing  effort  has  been  directed  to 
address  thermal-piezoelectric-mechanical  response  [6, 7,  8]. 

Many  engineering  smart  structures  have  one  dimension  much  smaller  than  the  other 
two  and  can  be  modeled  as  plates  if  there  are  no  initial  curvatures  associated  with  the  plane 
formed  by  the  two  large  dimensions.  The  present  capability  of  analyzing  thermopiezoelastic 
behavior  of  smart  plates  is  limited.  The  mathematical  models  are  generally  derived  from 
three-dimensional  (3-D)  elasticity  theory,  making  use  of  the  fact  that  the  thickness  is  small  in 
some  sense.  Most  analyses  prevailing  in  the  literature  are  so-called  ad  hoc  theories,  which 
can  be  generally  classified  as  Classical  Lamination  Theory  (CLT)  [9],  First-Order  Shear 
Deformation  Theory  (FOSDT)  [10],  Higher-Order  Theory  [8]  and  Layerwise  Theory  [11]. 
Layerwise  theory  can  produce  reasonable  results  at  the  cost  of  complex  models  and  expensive 
computation.  All  the  other  ad  hoc  approaches  are  doomed  to  fail,  especially  for  stress 
prediction  through  the  thickness,  even  for  shells  of  moderate  thickness.  The  reason  is  that 
these  theories  assume  the  displacements  to  be  functions,  while  in  fact  the  displacement 
field  of  a  layered  plate  may  have  discontinuous  derivatives  through  the  thickness. 

From  a  mathematical  point  of  view,  the  approximation  in  the  analysis  stems  from 
elimination  of  the  thickness  coordinate  from  the  independent  variables  of  the  governing 
equations  of  motion  for  the  plate  structure.  This  sort  of  approximation  is  inevitable  if  one 
wants  to  take  advantage  of  the  smallness  of  the  thickness  to  simplify  the  analysis.  However, 
other  approximations  that  are  not  absolutely  necessary  should  be  avoided.  For  example, 
for  small-strain  analysis  of  plates,  it  is  reasonable  to  assume  that  the  thickness,  h,  is  small 
compared  to  the  wavelength  of  deformation  of  the  reference  plane,  1.  However,  it  is  not  at  all 
reasonable  to  assume  a  priori  some  ad  hoc  displacement  field,  although  that  is  the  way  most 
existing  plate  theories  have  been  constructed. 

A  simple  and  accurate  model  of  composite  plates  and  shells,  namely.  Variational 
Asymptotic  Plate  and  Shell  Analysis  (VAPAS)  [12,  13,  14,  15],  was  developed  recently. 
VAPAS  starts  with  formulation  of  the  3-D  anisotropic  elasticity  problem  in  which  the 
deformation  of  the  reference  plane  is  expressed  in  terms  of  intrinsic  two-dimensional  (2-D) 
variables.  The  intrinsic  formulation  allows  the  body  to  undergo  arbitrarily  large  displacements 
and  global  rotations  subject  only  to  the  restriction  that  the  generalized  2-D  strains  are  small. 
The  Variational  Asymptotic  Method  (VAM)  [16]  is  then  used  to  systematically  reduce  the 
dimensionality  of  the  problem  by  taking  advantage  of  the  small  parameters  inherent  in  the 
problem  The  original  nonlinear  3-D  problem  is  thus  mathematically  split  into  a  linear  one¬ 
dimensional  (1-D)  through-the-thickness  analysis  and  a  nonlinear  2-D  plate/shell  analysis 
accounting  for  transverse  shear  deformation.  The  through-the-thickness  analysis  is  solved 
by  finite  element  method  and  provides  a  constitutive  model  between  the  generalized,  2-D 
strains  and  stress  resultants  as  well  as  recovery  relations  to  accurately  approximate  the  3-D 
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displacement,  strain  and  stress  fields  in  terms  of  2-D  variables  calculated  in  the  2-D  plate/shell 
analysis.  Numerical  examples  presented  in  previous  works  [12, 13,  .14, 15]  have  demonstrated 
for  mechanical  and  thermal  loading  that  although  the  resulting  theory  is  as  simple  as  a  single¬ 
layer  FOSDT,  the  recovered  3-D  displacement,  strain  and  stress  results  have  an  accuracy 
comparable  to  that  of  higher-order,  layer-wise  theories  with  many  more  degrees  of  freedom. 

The  present  work  extends  VAPAS  so  that  the  thermopiezoelastic  effects  of  smart  plates 
can  be  treated  in  the  same  framework.  Because  the  procedure  is  quite  similar,  the  authors 
have  chosen  to  repeat  some  formulae  and  text  from  their  previous  publications  in  order  to 
make  the  present  paper  more  self-contained.  The  present  theory  has  been  implemented  into 
the  computer  program  VAPAS  and  now  one  can  use  this  program  along  with  some  2-D 
plate  solver  (e.g.,  a  standard  plate  finite  element  in  commercial  finite  element  software  or 
in  a  flexible  multi-body  code  such  as  DYMORE  [17])  to  carry  out  an  accurate  and  efficient 
thermopiezoelastic  analysis  for  smart  composite  plates  made  with  piezoelectric  material. 

2.  3-D  Formulation 

A  point  in  the  plate  can  be  described  by  its  Cartesian  coordinates  Xi  (see  Figure  1),  where 
Xc,  are  two  orthogonal  lines  in  the  reference  plane  and  x^  is  the  normal  coordinate.  (Here 
and  throughout  the  paper,  Greek  indices  assume  values  1  and  2  while  Latin  indices  assume  I, 
2,  and  3.  Repeated  indices  are  summed  over  their  range  except  where  explicitly  indicated.) 
Letting  bj  denote  the  unit  vector  along  Xi  for  the  undeformed  plate,  one  can  then  describe  the 
position  of  any  material  point  in  the  undeformed  configuration  by  its  position  vector  f  from  a 
fixed  point  O,  such  that 

f(xi,a:2,a:3)  =  r(a;i,a;2)  (0 

where  r  is  the  position  vector  from  O  to  the  point  located  by  Xa  on  the  reference  plane. 
When  the  reference  plane  of  the  undeformed  plate  coincides  with  its  middle  plane,  it  naturally 
follows  that 

(f(a'i,a:2,a'’3))  =  r(xi,a:2)  (2) 

where  the  angle-brackets  denote  the  definite  integral  through  the  thickness  of  the  plate  and 
will  be  used  throughout  the  rest  of  the  development. 

When  the  plate  deforms,  the  particle  that  had  position  vector  r  in  the  undeformed  state 
now  has  position  vector  R  in  the  deformed  plate.  The  latter  can  be  uniquely  determined  by 
the  deformation  of  the  3-D  body.  Analogous  to  b,  for  the  undeformed  state,  another  triad  Bj 
is  introduced  for  the  deformed  configuration.  The  relation  between  Bj  and  bj  can  be  specified 
by  an  arbitrarily  large  rotation  specified  in  terms  of  the  matrix  of  direction  cosines  C{;X\,X2) 
so  that 

Bj  =  C'j,b,  Cij  =  Bihj  (3) 

subject  to  the  requirement  that  Bj  is  coincident  with  bj  when  the  structure  is  Undeformed. 
Now  the  position  vector  R  can  be  re  presented  as 

R(xi,X2,X3)  =  R(Xi,X2)  -|-X3B3(Xi,X2)  +  lUt(xi,X2,X3)Bj(xi,X2)  (4) 
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Figure  1.  Schematic  of  plate  deformation 

where  Wi  is  the  warping  of  the  normal-line  element.  In  the  present  work,  the  form  of  the 
warping  Wi  is  not  assumed,  as  in  most  plate  theories.  Rather,  these  quantities  are  treated 
as  unknown  3-D  functions  and  will  be  solved  for  later.  Equation  (4)  is  six  times  redundant 
because  of  the  way  warping  introduced.  Six  constraints  are  needed  to  rnake  the  formulation 
unique.  The  redundancy  can  be  removed  by  choosing  appropriate  definitions  of  R  and  Bj. 
One  can  define  R  similarly  as  Equation  (2)  to  be  tfte  average  position  through  the  thickness, 
from  which  it  follows  that  the  wafping  functions  must  satisfy  the  three  constraints 

{Wi{xi,X2,X3))  =0  (5) 

Another  two  constraints  can  be  specified  by  taking  B3  as  the  normal  to  the  reference  plane  of 
the  deformed  plate.  It  should  be  noted  that  this  choice  has  nothing  to  do  with  the  Kirchhoff 
hypothesis.  Indeed,  it  is  only  for  convenience  in  the  derivation.  In  the  Kirchhoff  assumption, 
no  local  deformation  of  the  transverse  normal  is  allowed.  However,  according  to  the  present 
scheme  we  allow  all  possible  deformation,  classifying  all  deformation  other  than  that  ot 
classical  lamination  theory  (CLT)  as  warping,  which  is  assumed  to  be  small.  This  assumption 
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is  valid  if  the  strain  is  small  and  if  the  local  rotation  {i.e.  the  rotation  of  the  normal  line  caused 
by  warping)  is  no  larger  than  the  order  of  the  strain  [18], 

Based  on  the  concept  of  decomposition  of  rotation  tensor  [19,  18],  the  Jauman-Biot- 
Cauchy  strain  components  for  small  local  rotation  are  given  by 

ry  =  i(Fi,  +  J=;,)-4  (6) 

where  is  the  mixed-basis  component  of  the  deformation  gradient  tensor  such  that 

Fij  =  BiGkg'^hj  (7) 


Here 


G/b  = 


dR 

dxk 


is  the  covariant  basis  vector  of  the  deformed  configuration  and  g^’  the  contravariant  base  vector 
of  the  undeformed  configuration  and  =  gA:  =  t>A  -  One  can  obtain  Ga:  with  the  help  of  the 
definition  of  so-called  generalized  2-D  strains  similarly  as  [20],  given  by 


R.a  —  Bq  -f  tQjjBj  (8) 

Bj.a  =  (— A'qjB^  X  B3  +  FasBa)  X  Bi  (9) 


where  and  Ka3  are  the  2-D  generalized  strains  and  comma  denotes  the  differentiation 
with  respect  to  the  coordinates.  . Here  one  is  free  to  set  £12  =  £21.  i  e. 

Bi-R2  =  B2-Ri  (10) 


which  can  serve  as  another  constraint  to  specify  the  deformed  configuration. 

With  the  assumption  that  the  strain  is  small  compared  to  unity,  which  has  the  effect  of 
removing  all  the  terms  that  are  products  of  the  warping  and  the  generalized  strains,  one  can 
express  the  3-D  strain  field  as 

r  =  F/jU;  +  FjC -h  F/j  -h  ri2'U)^2  (10 

where 


F  =  [Fn  2Fi2  F22  2Fi3  2F23  F33J’’ 

w  =  [rci  W2 

e  =  [£ii  2£i2  £22  Fii  Ki2  +  K2I  ^22]^ 
and  all  the  operators  are  defined  as: 


0 

0 

0 

’  1 

0 

o' 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

d 

dx:i 

0 

0 

0 

0 

1 

0 

d 

dx3 

0 

0 

0 

0 

0 

0 

d 

dx3  , 

0 

0 

0 
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In  the  present  work,  we  only  study  the  actuated  effect  caused  by  applied  thermal  or 
electrical  loads,  which  means  there  is  only  one-way  thermopiezoelastic  coupling.  The  changes 
of  temperature  and  electrical  fields  caused  by  deformation  of  the  plate  are  negligible,  and  the 
interactions  between  temperature  and  electric  field  are  not  considered  either.  Then  we  can  use 
the  Gibbs  free-energy  function  [21]  without  the  quadratic  terms  involving  temperature  and/or 
electric  field  to  carry  out  the  analysis.  The  free  energy  per  unit  area  (which  is  the  same  as  the 
tree  energy  of  the  normal-line  element)  can  be  written  as 

c.'  =  D  r  -  r^i:>Q  t  -  v'^Dd  (i?) 

where  T  is  the  difference  of  temperature  inside  the  structure  with  respect  to  the  reference 
temperature  when  the  plate  is  stress  free;  S  is  the  electric  field  vector;  D  is  the  3-D  6x6 
material  matrix,  which  consists  of  elements  of  the  elasticity  tensor  expressed  in  the  global 
coordinate  system  x,;  q  is  a  6x1  column  matrix  representing  the  3-D  thermal  expansion 
coefficients;  and  d  is  a  6x3  matrix  representing  the  3-D  strain-piezoelectric  coefficients. 
These  matrices  are  in  general  fully  populated.  However,  if  it  is  desired  to  model  laminated 
composite  plates  in  which  each  lamina  exhibits  a  monoclinic  symmetry  about  its  own  mid¬ 
plane  and  is  rotated  about  the  local  normal  to  be  a  layer  in  the  composite  laminated  plate, 
some  parts  of  D  will  always  vanish  [12]  no  matter  what  the  layup  angle  is,  and  a  and  d  will 
assume  the  following  forms: 

a  =  [ail  2ai2  a22  0  0  a33j^  (18) 

0  0  dii3 

0  0  2di23 

.  _  0  0  ^223 
2di3i  2di32  0 

2^231  2(^232  0 

0  0  0(333 

To  deal  with  applied  mechanical  loads,  we  will  at  first  leave  open  the  existence  of  a 
potential  energy  and  develop  instead  the  virtual  work  of  the  applied  mechanical  loads.  The 
virtual  displacement  is  taken  as  the  Lagrangean  variation  of  the  displacement  field,  such  that 

dR  =  6qgiBi  +  xsS^Bi^i  X  Ba  -I-  5wiBi  +  x  (20) 

where  the  virtual  displacement  of  the  reference  plane  is  given  by 

^gi=5u-Bf  (21) 

and  the  virtual  rotation  of  the  reference  plane  is  defined  such  that 

SBi  =  SipgjBj  X  Bi 


(22) 
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Since  the  strain  is  small,  one  may  safely  ignore  products  of  the  warping  and  the  loading  in  the 
virtual  rotation  term.  Then,  the  work  done  through  a  virtual  displacement  by  the  applied  loads 
r,Bi  at  the  top  surface  and  piBi  at  the  bottom  surface  and  by  the  body  force  fiBi  through  the 
thickness  is 


()TT  =  {Ti+^i  +  {4>i))SqQi  +  5ipQ^ 


2  {to,  —  Pa)  +  {^z4>o) 


+  (5  {Tiwf  +PiWi  +  (0iu;i))(23) 


where  r,.  Pi,  and  4>i  are  taken  to  be  independent  of  the  deformation,  ( )"*■  —  ( )  Ii3=h>  and 
()~  =  ( .  By  introducing  column  matrices  5q,  S'tp,  r,  P,  and  cf),  which  are  formed  by 
stacking  the  three  elements  associated  with  indexed  symbols  of  the  same  names,  and  using 
Equations  (1),  (3),  and  (4),  one  may  write  the  virtual  work  in  matrix  form,  so  that 

+  +  (24) 


SW  =  5qf  +  STl^m  +  S{i 


where 


/  =  T  +  ,5  +  (tA) 

■  h 
2 

m  =  I  j{t2-  P2)  +  (3:302) 


|(Ti  -  3l)  +  (x’30l) 


0 


(25) 


The  complete  statement  of  the  problem  can  now  be  presented  in  terms  of  the  principle  of 
virtual  work,  such  that 

(26) 


5U  -  SW  =  0 


In  spite  of  the  possibility  of  accounting  for  nonconservative  forces  in  the  above,  the  problem 
that  governs  the  warping  is  conservative  when  Ti,  pi,  and  0,-  are  taken  to  be  independent  of  the 
deformation.  Thus,  one  can  pose  the  problem  that  governs  the  warping  as  the  minimization 
of  a  total  potential  functional 

n  =  u  +  w  (27) 


so  that 

511  =  0  (28) 

in  which  only  the  warping  displacement  is  varied,  subject  to  the  constraints  Equation  (5). 
This  implies  that  the  potential  of  the  applied  mechanical  loads  for  the  functional  governing 
warping  is  given  by 

W  =  -  P'^w-  -  (0^iu)  (29) 

By  the  principle  of  minimum  total  potential  energy,  one  can  solve  for  the  unknown  warping 
functions  by  minimizing  the  functional  in  Equation  (27)  subject  to  the  constraints  of  Equation 
(5).  Up  to  this  point,  this  is  simply  an  alternative  formulation  of  the  original  3'D  elasticity 
problem.  If  we  attempt  to  solve  this  problem  directly,  we  will  meet  the  same  difficulty  as 
solving  any  full  3-D  elasticity  problem.  Fortunately,  the  VAM  can  be  used  to  calculate 
the  3-D  warping  functions  asymptotically.  Although  through-the-thickness  analysis  is  one 
dimensional  and  can  be  solved  analytically  [  12],  we  prefer  to  use  finite  element  discretization 
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to  solve  the  minimization  problem  for  the  sake  of  dealing  with  multiple  layers  and  arbitrary 
monoclinic  material,  and  connecting  with  any  2-D  standard  plate  solver  which  is  normally 
implemented  using  finite  element  method.  A  5-noded  isoparametric  element  is  used  because 
we  need  the  second-order  warping  functions  for  the  purpose  of  recovering  the  original  3-D 
fields  up  to  the  second  order  of  h/l.  These  second-order  warping  functions  are  piecewise, 
fourth-degree  polynomials  as  shown  in  [12],  Discretizing  the  transverse  normal  line  into  1-D 
finite  elements,  one  can  express  the  warping  field  as 

=  5(x3)V'(xi,a:2)  (30) 

where  S  is  the  shape  function  and  V  is  the  nodal  value  of  warping  field  along  the  transverse 
normal.  Substituting  Equation  (30)  into  Equation  (27),  one  can  express  the  total  energy  in 
discretized  form  as 

2n  =  V'^EV  +  -f-  Dm.Vx  +  Dhi.V^f) 

+  +  VjDi.i.V,  +  V^DumV^ 

■  +2{VlDi,,e  +  VlDi,,e  +  VjDi,i,V,2) 

-  21^^07,  -  2e^a,  -  2V[a;,  -  2VJq,2 

-  2V'^8h  -  2fSe  -  2y[4  -  2V;[4  -f-  2V'^L  (31) 

where  L  contains  the  load  related  terms  such  that 

L  =  -S^'^r-S-'^^-\S^6)  (32) 

The  new  matrix  variables  carry  the  properties  of  both  the  geometry  and  material: 


E;=([45FD[r,5]) 

Dm  =  {[ThSfDr,) 

Dm,  =  ([rA5]^D[r,,5]) 

Dm,  =  ([r,,5]^/?[45}) 

A.  =  {rjDT,) 

Dta,  =  ([r,,5]^/7[r,,5]) 

Di,i,  =  {{rt,sf  D[ri,s]) 

Di,i,  =  {nsf  D[rt,S]) 

Di,e  =  {[Vi^sfor,) 

Di,,  =  ([45]^Dr,) 

aH  =  {[rhSfDaT) 

a*  ={rjDaT) 

ai,  =  {[Ti,Sf  Dci  t) 

cci,  =  {[Fi.Sf  Da  T) 

SH^\[ThSYDdS)  S,  ={rjDd£) 

£k^{[ri,SfDd£)  4  =^{[ri,S]'^Dd£)  (33) 

Although  the  theory  itself  allows  for  an  exact  representation  for  arbitrary  temperature  and 
electric  field  distribution  through  the  thickness,  here  fourth-degree  polynomials  are  used 
to  approximate  both  distributions  for  each  normal-line  element.  The  discretized  form  of 
Equation  (5)  is 

V'^Hip  =  0  (34) 

where  H  =  and  ip  is  the  normalized  kernel  matrix  of  E  such  that  ip'^Hip  =  /.  Now 

our  problem  is  transformed  to  minimize  Equation  (31)  numerically,  subject  to  the  constraints 
in  Equation  (34). 
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3.  Dimensional  Reduction 

To  rigorously  reduce  the  original  3-D  problem  to  a  2-D  plate  problem,  one  must  attempt  to 
reproduce  the  energy  stored  in  the  3-D  structure  in  a  2-D  formulation.  This  dimensional 
reduction  can  only  be  done  approximately,  and  one  way  to  do  it  is  by  taking  advantage  of  the 
smallness  of  h/l.  The  small  parameter  e,  representing  the  order  of  the  generalized  2-D  strains 
e  has  already  been  taken  advantage  of  when  we  derive  Equation  ( 1 1 ).  To  reduce  the  number  of 
small  parameters  in  the  asymptotic  analysis,  it  is  reasonable  to  assume  that  the  strains  caused 
by  temperature  and  electricity  are  of  the  order  s.  Thus,  the  quantities  of  interest  assume  the 
following  orders; 

f3^p{h/iy£  fa  p{h/l)s 

nxa  ~  ph{h/l)£  OiT  ~  e  ~  f  •  (35) 

where  p  is  the  order  of  the  elastic  constants  (all  of  which  are  assumed  to  be  of  the  same  order). 

Having  assessed  the  orders  of  all  the  interested  quantities,  the  VAM  can  be  used  to 
mathematically  perform  a  dimensional  reduction  of  the  3-D  problem  to  a  series  of  2-D  models. 
This  method  requires  one  to  find  the  leading  terms  of  the  functional  according  to  the  different 
Orders.  Since  only  the  warping  is  varied,  the  leading  terms  needed  are  all  of  those  terms 
associated  with  warping.  For  the  zeroth-order  approximation,  these  leading  terms  of  Equation 
(31)  are 

2ns  =  +  2V'^D,„e  -  2V^ah  -  2V^Sh  (36) 

The  zeroth-order  warping  functions  which  minimize  the  above  functional  subject  to 
constraints  Equation  (34)  can  be  obtained  by  usual  procedure  of  calculus  of  variation  as: 

V  =  Vo€  +  Vt  +  V£:=Vo  (37) 

Substituting  Equation  (37)  back  into  Equation  (31),  one  can  obtain  the  total  energy 
asymptotically  correct  up  to  the  order  of  as 

2no  =  -  2e^AV  -  (38) 

with 

A  =  {V^Dh,  +  D,,) 

■  NT  =  a,  +  ^{Vo"aH-DlVT) 

Ne  =  £,  +  ^iV^SH-DlV£)  (39) 

Although  the  energy  of  this  approximation  coincides  with  the  classical  plate  theories  for 
thermopiezoelastic  analysis,  we  have  not  used  any  ad  hoc  kinematic  assumptions  such  as 
the  Kirchhoff  assumption  to  obtain  this  result.  Moreover,  the  transverse  normal  strain  from 
our  zeroth-order  approximation  is  not  zero. 

It  is  understood  that  our  zeroth-order  approximation  will  give  the  same  stress  results  as 
what  is  obtained  from  CLT,  i.e.,  all  the  transverse  stress  components  which  are  very  important 
for  analyzing  the  failure  of  composite  plates  cannot  be  predicted.  One  must  carry  out  the 
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next  approximation  so  that  those  quantities  can  be  approximately  predicted.  To  obtain  the 
first-order  approximation,  we  simply  perturb  the  zeroth-order  result,  such  that 

V  =  Vo  +  Vi  (40) 

Substituting  Equation  (40)  back  into  Equation  (1 1)  and  then  into  Equation  (31),  one  can  obtain 
the  leading  terms  for  the  first-order  approximation  as 

2nj  —  EVi  -f-  2V^ -h  TV2  -^2^,2  +  2V^ Lx  +  2V^ -F  2V^ L  (41) 

with 

Lx  =  (D/i/,  —  DJiJVx,i  +  {Dhi2  —  +  ocii,i  +  0:12,2 

~  i  -}-  {Dhl2  —  D1i^)Vs,2  +  -I-  5/2,2 

D,  =  {Dm,  -  DljVo  -  D,,, 

Do  =  {Dm2  -  Dl^)Vo  -  D,,,  (42) 

Integration  by  parts  with  respect  to  the  in-plane  coordinates  is  used  here  and  hereafter 
whenever  it  is  convenient  for  the  derivation,  because  the  present  goal  is  to  obtain  an  interior 
solution  for  the  plate  without  consideration  of  edge  effects.  Note  that  the  treatment  of  edge 
effects  is  itself  a  very  important  problem  to  tackle,  but  it  is  outside  the  scope  of  the  present 
work. 

Similarly  as  in  the  zeroth-order  approximation,  one  can  solve  the  first-order  warping  field 
as 

=  Vjie.i  +  V'i2f  2  +  Vix  +  Vis  +  V'li  (43) 

and  obtain  a  total  energy  that  is  asymptotically  correct  up  to  the  order  of  p{h/iys,  given  by 
2ni  =  e'^Ae  +  e^5e,i  -f  2eJCe,2  +  2  -  2e^F  -  2e^Fx  -  2e^ Fe  (44) 

where 

B  =  VfA,„Vb  +  ViTA 
C  =  ViD,,,,V„  +  \(\f,D.,  +  DjVn) 

'P  =  VoD,,,,Vo  +  V,ID, 

F  =  L  —  +  VuF,i  +  D2V1L2  +  ^12^,1)  (45) 

with  the  non-mechanical  load  due  to  temperature 

Fx  =  Nx  +  Vq  Di,i,Vx,n  +  Vq  A2/2^,22  +  Vq  {Dip^  +  Df^jjVx^u  (46) 
+  +  ^\2^t,2  +  D^Vix,\  +  DIVix,2) 

and  the  non-mechanical  load  due  to  electric  field 

Fe  =  Ne-  Vq  Dhi^Ve^y  -  Vq  Dhi2Vs,2  -  D^i^Ve^i  -  D^i^Ve^2 
+  VcTAi/iVf.n  +  V^Di2i2Vs,22  +  V^  {Di,i2  +  Df^i^)V£^i2 
+  +  yi2^e,2  +  D1Vis,i  +  T)^Vi£:,2) 


(47) 
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Here  the  monoclinic  symmetry  has  already  been  used  to  obtain  the  asymptotically  correct 
energy  in  Equation  (44).  The  applied  mechanical  loads,  temperature  and  electric  field  should 
not  vary  rapidly  over  the  plate  surface,  so  that  F,  Ft  and  Fs  will  be  of  sufficiently  high  order 
to  meet  the  requirement  of  asymptotical  correctness. 

4.  Transforming  into  a  Reissner-Mindlin  Model 

Although  Equation  (44)  is  asymptotically  correct  through  the  second  order  and  use  of  this 
free  energy  expression  is  possible,  it  involves  more  complicated  boundary  conditions  than 
necessary  since  it  contains  derivatives  of  the  generalized  strain  measures.  To  obtain  an  energy 
functional  that  is  of  practical  use,  one  can  transform  the  present  approximation  into  a  Reissner- 
Mindlin  type  model.  It  should  be  noted  that  fitting  the  asymptotic  energy  into  such  model 
is  just  one  possible  choice,  and  the  possibility  of  fitting  the  same  energy  into  other  more 
sophisticated  2-D  plate  models  exists. 

In  a  Reissner-Mindlin  model,  there  are  two  additional  degrees  of  freedom,  which  are 
the  transverse  sheaf  strains.  These  are  incorporated  into  the  rotation  of  transverse  normal.  If 
we  introduce  another  triad  B*  for  the  deformed  Reissner-Mindlin  plate,  the  definition  of  2-D 
strains  becomes 

B-.  =  (-/CjB;  X  B;  +  KWBl)  x  B*  (48) 

where  the  transverse  shear  strains  are  7  =  [2713  2723]^.  One  can  express  the  classical  strain 
measures  e  in  terms  of  the  strain  measures  of  the  Reissner-Mindlin  plate  model  as 

e  =  n-  Da7,o  (49) 

where 

r  0  0  0  1  0  0  1^ 

^“[000010 

^  [  0  0  0  0  1  0  ]^ 

^[  0  0  0  0  0  1 

%  =  [e;,  4,  /(•;,  iCIj+A-J,  (50) 

Now  one  can  express  the  energy.  Equation  (44),  correct  to  second  order,  in  terms  of  strains  of 
the  Reissner-Mindlin  model  as 

2ni  =  'RFATI  -  27^^ADa7,a  +  +  27^^C7^,2  +  7^5^)7^,2 

-  2V7F  -  2R^Ft  +  2')'lDaNT 

-2R^Fs  +  2-i^,D^Ne  (51) 

The  generalized  Reissner-Mindlin  model  used  in  the  2-D  thermopiezoelastic  analysis  is  of  the 
form 


2117?,  =  A'R.  —  2TZ^{FTi  +  FTR.  +  P'£Tz)+')^G'Y  —  2'y'^{F.y  +  FT'y  +  F£y)(52) 
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To  find  an  equivalent  Reissner-Mindlin  model  Equation  (52)  for  Equation  (51),  one  has 
to  eliminate  all  partial  derivatives  of  the  classical  2-D  strain  measures.  The  equilibrium 
equations  are  used  to  achieve  this  purpose.  From  the  two  equilibrium  equations  balancing 
bending  moments  with  applied  moments  nia  which  is  calculated  from  Equation  (25),  one  can 
obtain  the  following  formula 

V^iAnc  -  Fn,a  -  FTn,a  -  Fen,a)  =  G'y-F^-  Ft,  "  -  |^^  }  (53) 

Using  Equation  (53),  one  can  rewrite  Equation  (51)  as 

2ni  =  ATlW  01-211^  {F^FT^Fe)-2^'^V^NT,a-2-fV^Ne,^+U\5A) 

where 

U*  =  TllBUi  +  27^^C7^,2  +  (55) 

and 

B  =  B  +  AV^G-'^V'l  A 
C  =  C  +  AViG-'^VlA 

D  =  D  +  AV2G-^VIA  (56) 

If  we  can  drive  U*  to  zero  for  any  IZ,  then  we  have  found  an  asymptotically  correct  Reissner- 
Mindlin  plate  model.  For  generally  anisotropic  plates,  however,  this  term  cannot  be  driven 
to  zero;  but  we  can  minimize  the  error  to  obtain  a  Reissner-Mindlin  model  that  is  as  close  to 
asymptotical  correctness  as  possible.  The  accuracy  of  the  Reissner-Mindlin  model  depends 
on  how  close  to  zero  one  can  drive  this  term  of  the  energy. 

One  could  proceed  with  the  optimization  at  this  point,  but  the  problem  will  require  a 
least  squares  solution  for  3  unknowns  (elements  of  the  shear  stiffness  matrix  G)  from  a 
linear  system  of  78  equations  (12x12  and  symmetric),  a  very  rigid  optimization  problem. 
The  solution  will  be  better  if  we  can  bring  more  unknowns  into  the  problem.  There  is  no 
unique  plate  theory  of  a  given  order  [22].  One  can  relax  the  constraints  in  Equation  (5)  to 
be  {wi)  =  constant  and  still  obtain  an  asymptotically  correct  strain  energy.  Since  the  zeroth- 
order  approximation  gives  us  an  asymptotic  model  corresponding  to  classical  plate  theory,  we 
only  relax  the  constraints  for  the  first-order  approximation.  This  relaxation  will  modify  the 
warping  field  to  be 

y  1  ==  Uiie,i  -+-  Viae, 2  +  Vil  +  ViT  +  Vu  +  +  L2e,2  (57) 

where  Li ,  L2  consist  of  24  constants.  The  remaining  energy  U*  will  also  be  modified  to  be 
u*  =  n^Bni  +  27^^C7^,2  +  nlbuo  (58) 

and 

B=B  +  2LjDi 
C=C+iL^D2+Di^L2) 
i)  =D  -(-  2L2D2 


(59) 
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Since  now  we  have  27  unknowns,  the  optimization  is  much  more  fljexible.  It  can  give  us  a  more 
optima]  solution  for  the  shear  stiffness  matrix  G  to  fit  the  second-order,  asymptotically-correct 
energy  into  a  Reissner-Mindlin  model.  In  other  words,  when  one  carries  out  the  optimization 
as  described  here,  one  finds  the  Reissner-Mindlin  model  that  describes  as  closely  as  possible 
the  2-D  energy  that  is  asymptotically  correct  through  the  second  order  in  h/l.  However,  the 
asymptotical  correctness  of  the  warping  field  to  that  same  order  can  only  be  ascertained  after 
obtaining  another  higher-order  approximation,  which  will  be  discussed  in  the  next  section. 

After  minimizing  U*,  the  “best”  total  energy  to  be  used  for  the  2-D  plate  Reissner- 
Mindlin  model  for  the  purposes  of  thermopiezoelastic  analysis  can  be  expressed  as: 

2Un  =  n^ATZ  -  2n'^{F  +  Ft  -h  Fs)  -t-  y^G'y  -  27^Pa(A'r,a  +  Ns, a)  (60) 

5.  Recovering  Relations 

From  the  above,  we  have  obtained  a  Reissner-Mindlin  plate  model  which  is  as  close  as 
possible  to  being  asymptotically  correct  in  the  sense  of  matching  the  total  energy.  The  stiffness 
matrices  A,  G,  load-related  terms  and  non-mechanical  stress  resultants  can  be  used  as  input 
for  a  plate  theory  derived  from  the  total  energy  obtained  here.  The  nonlinear  theory  presented 
in  [20]  is  an  appropriate  choice,  but  some  modification  of  the  loading  terms  is  required. 

However,  while  it  is  necessary  to  accurately  calculate  the  2-D  displacement  field  of  the 
plates,  this  is  not  sufficient  in  many  applications.  Ultimately,  the  fidelity  of  a  reduced-order 
model  such  as  this  depends  on  how  well  it  can  predict  the  3-D  results  in  the  original  3- 
D  problem.  Hence,  recovering  relations  should  be  provided  to  complete  the  reduced-order 
model.  By  recovering  relations,  we  mean  expressions  for  3-D  displacement,  strain,  and  stress 
fields  in  terms  of  2-D  quantities  and  x^.  For  validation,  results  obtained  for  the  3-D  field 
variables  from  the  reduced-order  model  must  be  compared  with  those  of  the  original  3-D 
model. 

For  an  energy  that  is  asymptotically  correct  through  the  second  order,  we  can  recover 
the  3-D  displacement,  strain  and  stress  fields  only  through  the  first  order  in  a  strict  sense 
of  asymptotical  correctness.  Using  Equations  (1),  (3),  and  (4),  one  can  recover  the  3-D 
displacement  field  through  the  first  order  as 

[  C31  ' 

Uzd  — “^2(13- Xz\  C32  > -t- -F  S'Fi  (61) 

[  Czz  ~  1  ^ 

where  Uzd  is  the  column  matrix  of  3-D  displacements  and  u^d  is  the  plate  displacements.  Cy 
are  the  components  of  global  rotation  tensor  from  Equation  (3).  From  Equation  (11),  one  can 
recover  the  3-D  strain  field  through  the  first  order  as 

r  =  rhS{Vo-{-Vi)  +  T,e  +  T,,SVo,i3-ri,SVo,2  (62) 

Then,  one  can  use  the  3-D  constitutive  relation 
<j  =  DT  -DaT-DdS 
to  obtain  3-D  stresses  <Ty . 


(63) 
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Since  we  have  obtained  an  optimum  shear  stiffness  matrix  G,  some  of  the  recovered 
3-D  results  through  the  first  order  are  better  than  classical  theory  and  conventional  FOSDT. 
(Note  that  conventional  FOSDT  has  no  rational  way  to  find  the  shear  stiffness  coefficients.) 
However,  for  the  transverse  normal  component  of  strain  and  stress  (Le.  Fas  and  <733),  it  can  be 
shown  that  the  agreement  is  not  satisfactory.  Let  us  recall  that  the  Reissner-Mindlin  theory  that 
has  been  constructed  only  ensures  a  good  fit  with  the  asymptotically  correct  3-D  displacement 
field  of  the  first  order  (while  energy  is  approximated  to  the  second  order).  Thus,  in  order  to 
obtain  recovering  relations  that  are  valid  to  the  same  order  as  the  energy,  the  VAM  iteration 
needs  to  be  applied  one  more  time. 

Using  the  same  procedure  listed  in  previous  section,  the  second-order  warping  can  be 
obtained  and  expressed  symbolically  as 


•  V2  =  V2ie,n  +  V22e, 12  +  ^3^,22  +  F2L  +  Far  +  F2£  (64) 

Equation  (64)  is  obtained  by  taking  the  original  first-order  warping  Vi  to  be  the  result  of  the 
first-order  approximation.  It  is  clear  that  V2  is  one  order  higher  than  Vi  which  confirms  that  Vi 
is  the  first-order  approximation.  One  might  be  tempted  to  use  Vi  in  the  recovering  relations. 
However,  the  VAM  has  split  the  original  3-D  problem  into  two  sets  of  problems.  As  far  as 
recovering  relations  concerned,  it  is  observed  that  the  normal-line  analysis  can  at  best  give  us 
an  approximate  shape  of  the  distribution  of  3-D  results.  The  2-D  plate  analysis  will  govern 
the  global  behavior  of  the  stracture.  Since  V 1  is  the  warping  that  yields  a  Reissner-Mindlin 
plate  model  that  is  as  close  as  possible  to  being  asymptotically  correct,  we  should  still  use  Vi 
in  the  recovering  relations  instead  of  Iq.  By  doing  this,  we  choose  to  be  more  consistent  with 
Reissner-Mindlin  plate  model  while  compromising  somewhat  on  the  asymptotical  correctness 
of  the  shape  of  the  distribution.  It  has  been  verified  by  numerical  examples  that  this  is  a  good 
choice. 

Hence,  we  write  the  3-D  recovering  relations  for  displacement  through  the  second  order 


as 


Usd  —  +  3:^3  \ 


^31 


)  -f  S{\'o  +  1^1  +  V2) 


(65) 


^32  : 

C33  —  1 

and  the  strain  field  through  the  second  order  is 

r  =  ThSiVo  -f  Fi  -h  14)  +  r,€  +  r,, 5(14.1  +  Fi,i)  q-  r,, 5(14,2 + Fi,2)  m 


Again  the  stresses  through  the  second  order  can  be  obtained  from  the  3-D  constitutive  law. 
Equation  (63). 


6.  Numerical  Examples 

The  computer  program  VAPAS  has  been  extended  to  implement  the  present  theory.  Several 
numerical  examples  are  given  here  to  validate  the  proposed  theory  and  code  against  the  3-D 
thermopiezoelasticity  solutions  with  one-way  coupling  that  are  specialized  from  [23]. 

First  to  assess  the  asymptotical  correctness  of  the  proposed  theory,  we  study  a  cylindrical 
bending  type  problem  for  a  single-layer  plate  made  with  a  piezoelectric  material  with  isotropic 
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mechanical  properties  (E  as  the  Young’s  modulus,  u  Poisson’s  ratio)  and  rf  as  the  strain- 
piezoelectric  constants  in  the  both  directions  of  the  reference  plane.  The  plate  is  simply 
supported  with  width  L  along  xi  axis  (the  “lateral”  direction)  and  infinitely  long  along 
the  .r2  axis  (the  “longitudinal”  direction)  under  the  following  electricity  applied  in  voltage 
throughout  the  plate: 

<p  =  4>osm(J^^  (67) 

Let  us  assume  the  thickness  of  the  plate  is  h,  and  the  normalized  thickness  coordinate 
C  =  X3//1.  Then,  the  small  parameter  used  in  our  theory  is 

<5  =  -  =  —  (68) 

I  L 

Table  1  lists  the  nontrivial  displacements  and  stresses.  The  exact  solutions  are  obtained 
based  on  [23]  and  expanded  into  a  series  in  terms  of  5  with  o(*)  denoting  terms  asymptotically 
smaller  than  the  order  of  *.  The  present  theory  can  predict  the  correct  results  up  to  the  second 
order  of  6  with  respect  to  the  leading  terms  for  each  3-D  quantities,  which  clearly  demonstrate 
that  our  theory  is  asymptotically  correct  up  to  the  second  order  for  this  particular  problem, 
although  the  prediction  of  transverse  components  for  this  problem  is  out  of  the  power  of  our 
theory. 

However,  this  should  riot  mislead  the  reader  to  assume  that  the  present  theory  is 
asymptotically  correct  up  to  the  second  order  in  general.  The  authors  are  aware  that  the 
proposed  theory  can  be  at  best  asymptotically  correct  up  to  the  second  order  for  particular 
cases.  For  the  general  case,  however,  the  theory  can  only  be  interpreted  as  that  Reissner- 
Mindlin  model  which  is  the  closest  to  being  asymptotically  correct.  To  illustrate  the  above 
statement,  we  provide  the  results  for  the  same  piezoelectric  plate  under  a  transverse  surface 
mechanical  load  in  addition  to  the  aforementioned  electrical  charge: 


Table  2  only  lists  the  nontrivial  displacement  results.  One  can  observe  from  Table  2  that  there 
is  a  slight  difference  for  the  second  order  prediction  (relative  to  the  dominant  terms)  between 
the  present  theory  and  exact  solution.  It  is  interesting  to  note  that  if  one  sets  u  equal  to  zero 
the  difference  disappears.  Evidently  some  information  belonging  to  second  order  and  indeed 
included  in  the  asymptotically  correct  energy  cannot  be  captured  in  a  Reissner-Mindlin  type 
model.  When  we  transform  the  asymptotically  correct  model  into  a  Reissner-Mindlin  model, 
this  information  is  lost. 

However,  since  the  loss  is  small  in  comparison  to  the  dominant  terms,  the  numerical 
difference  between  the  present  theory  and  the  exact  solution  is  expected  to  be  small.  To  verify 
this  expectation,  we  will  present  some  numerical  results  for  piezoelectric  plates  to  demonstrate 
the  accuracy  of  our  theory.  We  study  a  single-layer  plate  with  h=l  mm  and  L=4  mm  under 
the  applied  electricity  as  in  Equation  (67)  with  (^o=100  V  and  mechanical  load  on  the  surfaces 
as  in  Equation  (69)  with  po=l  MPa.  The  piezoelectric  material  properties  are  taken  from  [8] 

=  Ex  —  63  GPa 
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Table  1.  Three-dimensional  displacement  and  stresses  under  uniform  temperature  change 
through  the  thickness 


Normalized  lateral  displacement 


1/+1)  c^i  7r(12C^-l)i/ 


7T  ^  24  '  5760(1/- 1) 

Present  +  o(S^) 


Normalized  transverse  displacement  i]^) 

Exact  K  +  +  o(^^) 

Present  uC  +  +  o(5^) 


Normalized  lateral  in-plane  stress 


Exact 

Present  +  o{S^) 


Normalized  longitudinal  in-plane  stress  (;^-) 


Exact  -1  -  +  o(S^) 

Present  -1  -  +  o(fi) 


Normalized  lateral  transverse  shear  stress  i-§^) 
Present  Not  available 


Normalized  lateral  transverse  shear  stress  (;^-) 


Ettact  +  o{S^) 

Present  Not  available 


Table  2.  Three-dimensional  displacements  under  linearly  distributed  temperature  change 
through  the  thickness 


Normalized  lateral  displacement  (Ui) 


12hpoC(i/^—l) 

Ett^ 


h{i^-\-l)\poC\{20iy-^0)C+9u^e)-10Ed<Po\ 


Present 

Normalized  transverse  displacement  (C/s) 


Exact  tttt.M-=-t)j-4  _ 


Present  - 
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Figure  2.  Distribution  of  the  3-D  stress  gu  vs  the  thickness  coordinate.  Solid  line:  exact 
solution:  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Normalized  thickness  coordinate 


Figure  3.  Distribution  of  the  3-D  stress  (722  vs  the  tliickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 


Glt  ~  Gtt  “  24.6  GPa 

^LT  =  ^TT  0.28 

rfll3  =  ^223  =  150  X  10-^2  ^  (70) 

Figures  2-5  plot  the  nontrivial  components  of  3-D  stress  distribution  through  the  thickness. 
(Note  that,  because  the  2-D  variables  are  either  sine  or  cosine  functions  of  xi,  and  (T33  are 
plotted  for  the  position  xi  =  L/2,  and  (Tqs  are  plotted  for  the  position  xi  =  0  or  a’l  =  L.) 
One  can  observe  that  VAPAS  results  are,  almost  on  the  top  of  exact  solutions  and  much  better 
than  the  results  of  CLT  and  FOSDT.  The  loss  of  information  is  almost  negligible. 

The  present  theory  is  formulated  with  sufficient  generality  to  carry  out  a  thermopiezoe¬ 
lastic  analysis  for  arbitrary  composite  laminated  piezoelectric  plates  made  with  monoclinic 
material  with  a  computational  cost  equivalent  to  that  of  FOSDT.  To  demonstrate  this  fact,  we 
study  a  more  challenging  and  realistic  problem.  It  is  a  four-layer  smart  plate  with  h=\  mm 
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Figure  4.  Distribution  of  the  3-D  stress  <7iz  vs  the  thickness  coordinate.  Solid  line:  exact 
solution:  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Normalized  thickness  coordinate 


Figure  5.  Distribution  of  the  3-D  stress  azz  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 


and  Z/=10  mm  (see  Figure  6).  The  two  face  sheets  are  made  with  piezoelectric  material  with 
properties  as  in  Equation  (70)  and  the  inside  layers  are  normal  graphite/epoxy  composites 
with  the  following  properties: 

£;L  =  l72GPa  Et  =6.9GPa 
Glt  ~  3-4  GPa  Gtt  =  1'4  GPa 

t/xr  ~  0.25  utj-  =  0.25  (71) 

The  piezoelectric  layers  are  each  0.1  mm  thick,  and  the  regular  composite  layers  are  each 
0.4  mm  thick.  The  layup  scheme  is  [0°/  -45°/45°/0°]  from  bottom  to  top.  An  electric  charge 
according  to  Equation  (67)  with  (po=lO  V  is  applied  to  both  piezoelectric  layers  with  the 
positive  direction  align  with  the  xz  coordinate.  The  recovered  stresses  are  plotted  in  Figures 
7  “  12.  Again,  there  is  an  excellent  agreement  between  VAPAS  results  and  those  of  the  exact 
solution  except  the  transverse  normal  stress.  VAPAS  cannot  predict  a  very  accurate  transverse 


Simple  Thermopiezoelastic  Model  for  Composite  Plates 


19 


normal  stress  for  this  case  because  this  quantity  contains  significant  terms  with  order  higher 
than  the  second  order  which  is  beyond  the  capability  of  VAPAS,  as  we  observed  from  Table  1 . 
Nevertheless  VAPAS  qualitatively  captures  the  shape  of  the  distribution  through  the  thickness. 
Careful  readers  may  also  find  that  there  is  an  obvious  discontinuity  of  the  distribution  of  the 
transverse  normal  stress  which  should  not  be  there.  The  reason  for  this  violation  is  that  our 
theory  is  displacement-based  and  at  most  asymptotically  correct  up  to  the  second  order.  The 
approximation  causes  this  discontinuity.  However,  it  is  comforting  to  notice  that  the  actual 
values  for  this  stress  component  are  very  small  relative  to  the  rest. 

To  investigate  the  behavior  of  a  piezoelectric  plate  under  a  combination  of  different  kinds 
of  loads,  we  apply  a  mechanical  load  according  to  Equation  (69)  with  po  =  IMPa  onto  the 
surfaces.  All  the  components  of  the  stress  due  to  the  combination  of  these  two  loads  are 
plotted  in  Figures  13  -  18.  This  time  VAPAS  achieves  excellent  agreement  with  the  exact 
solution  for  all  six  stress  components,  including  the  transverse  normal  stress.  This  happens 
because  stresses  caused  by  mechanical  loads  are  of  orders  lower  than  the  second  order  and 
within  the  power  of  the  present  theory. 

The  thermoelastic  behavior  of  composite  plates  has  been  studied  in  [15].  The  current 
version  of  VAPAS  can  reproduce  all  the  results  there.  No  additional  examples  will  be  given 
here  for  the  sake  of  brevity. 
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Figure  7.  Distribution  of  the  3-D  stress  cth  vs  the  thickness  coordinate.  Solid  line:  exact 
solution:  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Figure  8.  Distribution  of  the  3-D  stress  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dasli/short-dash  line:  CLT. 


Figure  9.  Distribution  of  the  3-D  stress  &22  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Normalized  thicknesscoordinate 

Figure  10.  Distribution  of  the  3-D  stress  crjs  vs  the  thickness  coordinate.  Solid  line:  exact 
solution:  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 


Figure  11.  Distribution  of  the  3-D  stress  C723  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/shorl-dash  line:  CLT. 


Figure  12.  Distribution  of  the  3-D  stress  ^33  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Figure  13.  Distribution  of  the  3-D  stress  crii  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Figure  14.  Distribution  of  the  3-D  stress  a\2  ys  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Figure  15.  Distribution  of  the  3-D  stress  0-22  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Figure  16.  Distribution  of  the  3-D  stress  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short- dash  line:  CLT. 


Figure  17.  Distribution  of  the  3-D  stress  a23  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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Figure  18.  Distribution  of  the  3-D  stress  (T33  vs  the  thickness  coordinate.  Solid  line:  exact 
solution;  dots:  VAPAS;  dashed  line:  FOSDT;  long-dash/short-dash  line:  CLT. 
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A  tliermopiezoelastic  model  of  composite  piezoeletric  plates  has  been  developed  by  use  of  the 
variational  asymptotic  method.  A  general  2-D  constitutive  law  that  is  as  close  to  asymptotical 
coiTectness  as  possible  has  been  obtained  by  solving  the  1-D  through-the-thickness  analysis. 
The  original  3-D  results  have  been  reproduced  as  accurately  as  possible  from  a  Reissner- 
Mindlin  type  plate  analysis.  Numerical  examples  show  that,  although  the  resulting  theory  is  as 
simple  and  efficient  as  a  single-layer  first  order  shear  deformation  theory,  it  can  approximate 
the  3-D  exact  solution  very  accurately. 

The  present  work  differs  from  previous  work  on  this  topic  in  the  literature  in  the 
following  five  aspects: 

(i)  The  present  formulation  is  in  an  intrinsic  form  which  is  suitable  for  both  geometrically 
nonlinear  and  linear  plate  theories. 

(ii)  Without  using  any  kinematical  ad  hoc  assumptions,  the  dimensional  reduction  from  3-D 
to  2-D  is  carried  out  systematically  by  using  the  variational  asymptotic  method. 

(iii) .  The  Reissner-Mindlin  type  model  proposed  in  the  present  theory  is  not  a  FOSDT. 

The  present  theory  differs  from  FOSDT  by  representing  all  the  deformations  that  are 
purposely  eliminated  in  the  development  of  FOSDT.  This  is  accomplished  through 
allowing  all  possible  deformation  in  the  3-D  warping  functions,  which  are  solved  in 
turn  by  the  variational  asymptotic  method. 

(iv)  All  the  loads  (mechanical,  thermal,  and  electric)  that  can  be  applied  to  a  piezoelectric 
plate  can  be  handled  by  the  present  theory.  Moreover,  all  the  load  distribution  through 
the  thickness  can  be  arbitrary  and  is  approximated  in  VAPAS  by  a  layerwise  fourth- 
degree  polynomial.  This  is  more  realistic  and  accurate  than  most  existing  models,  in 
which  only  a  constant  or  linear  distribution  is  allowed. 

(v)  The  present  theory  decouples  the  modeling  process  of  the  plate  completely  from  the  2- 
D  plate  problem  described  on  the  reference  plane  so  that  the  obtained  2-D  generalized 
constitutive  model  can  be  used  as  input  for  any  other  2-D  standard  plate  solver.  The 
flexibility  of  VAPAS  connecting  to  standard  plate  and  shell  solvers  can  help  the  structural 
analysts  focus  more  on  solving  2-D  problem  for  different  situations. 

The  computer  program  VAPAS  can  now  be  used  along  with  a  2-D  Reissner-Mindlin  type 
plate  solver  to  perform  an  efficient  yet  accurate  and  detailed  analysis  for  thermopiezoelastic 
behavior  of  laminated  piezoelectric  composite  plates.  Such  a  tool  will  be  very  useful  for 
designers  of  piezoelectric  plates  to  more  efficiently  carry  out  accurate  tradeoff  studies.  It  has 
to  be  emphasized  that  the  present  work  only  deals  with  thermal  and/or  piezoelectric  actuation 
of  smart  plates.  Changes  of  temperature  and  electric  fields  caused  by  deformation  of  the  plate, 
which  is  the  so-called  sensing  capability  of  smart  plates,  cannot  be  treated  by  means  of  the 
present  theory.  To  extend  the  theory  to  deal  with  this  class  of  problems  is  possible  but  requires 
significant  effort. 
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Abstract 

A  novel  integration  scheme  for  nonlinear  dynamics  of  geometrically  exact  shells  is  developed  based  on 
the  inextensible  director  assumption.  The  new  algorithm  is  designed  so  as  to  imply  the  strict  decay  of  the 
system  total  mechanical  energy  at  each  time  step,  and  consequently  unconditional  stability  is  achieved 
in  the  nonlinear  regime.  Furthermore,  the  scheme  features  tunable  high  frequency  numerical  damping 
and  it  is  therefore  stiffly  accurate.  The  method  is  tested  for  a  finite  element  spatial  formulation  of  shells 
based  on  mixed  interpolations  of  strain  tensorial  components  and  on  a  two-parameter  representation  of 
director  rotations.  The  robustness  of  the  scheme  is  illustrated  with  the  help  of  numerical  examples. 


1  Introduction 

The  formulation  of  integration  algorithms  for  nonlinear  dynamics  of  geometrically  exact  shells  is  the  focus 
of  this  work.  The  partial  differential  equations  governing  this  class  of  problems  are  known  to  present  a 
rich  mathematical  structure.  In  particular,  the  resulting  models  are  Hamiltonian  systems  characterized  by 
a  symplectic  nature  and  associated  with  conservation  laws  that  stem  from  symmetries  of  the  Hamiltonian, 
The  linear  and  angular  momentum  as  well  as  the  total  mechanical  energy  are  conserved  for  free  motions  of 
such  systems. 

The  understanding  of  the  geometric  characteristics  of  the  governing  equations  has  been  historically  con¬ 
fined  to  the  fields  of  anal5d:ical  mechanics  and  pure  mathematics.  Surprisingly,  this  knowledge  has  been 
seldom  used  for  the  development  of  numerical  methods.  Indeed,  the  study  of  new  integration  algorithms  has 
been  traditionally  preoccupied  with  the  development  of  methods  applicable  to  vast  classes  of  problems,  for 
example  the  class  of  differential/algebraic  equations,  or  hyperbolic  conservation  laws.  Consequently,  classical 
methods  rarely  preserve  the  underlying  structure  of  the  problem  being  solved,  and  hence,  such  structure  is 
lost  in  the  numerical  solution. 

This  approach  also  limits  the  possible  theoretical  analyses  of  the  schemes,  which  are,  more  often  than 
not,  confined  to  linear  or  model  cases.  For  instance,  it  is  customary  to  characterize  integration  schemes 
for  structural  dynamics  by  studying  their  behavior  when  applied  to  a  linear  oscillator.  This  approach  is 
clearly  not  adequate  when  dealing  with  highly  nonlinear  problems  such  as  the  dynamics  of  geometrically 
exact  shells. 

A  new  approach  to  the  design  of  integration  algorithms  attempts  to  bridge  the  divide  between  theoretical 
and  numerical  mechanics.  Under  this  new  paradigm,  numerical  schemes  are  “backward-engineered”  to 
preserve  some  important  qualitative  features  of  the  governing  equations.  Fittingly,  this  approach  is  now 
called  geometric  integration  in  the  mathematical  community  [20].  Attempts  at  designing  “geometry-aware” 
algorithms  for  structural  dynamics  problems  can  be  traced  back  to  the  work  of  Simo  and  co-workers  who 
analyzed  the  problems  of  rigid  body  dynamics  [32],  nonlinear  elasto-dynamics  [29],  geometrically  exact 
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shells  [30],  and  geometrically  exact  beams  [31].  In  all  cases,  the  idea  was  to  design  algorithms  that  ensure 
the  discrete  preservation  of  the  total  mechanical  energy  and  of  the  linear  and  angular  momenta  of  the  system. 

When  integrating  linear  and  nonlinear  finite  element  models,  the  implications  of  the  discrete  equations 
stiffness  must  be  carefully  considered.  Indeed,  high  frequencies  are  an  artifact  of  the  spatial  discretization 
process  and  do  not  reflect  the  high  frequencies  of  the  original  infinite  dimensional  problem.  The  need  for 
high  frequency  numerical  dissipation  has  been  recognized  in  the  past  for  linear  problems  [23].  When  dealing 
with  complex  nonlinear  systems,  numerical  dissipation  becomes  indispensable.  Indeed,  nonlinearities  provide 
a  mechanism  for  transferring  energy  from  low  to  high  frequency  modes.  Consequently,  numerical  solutions 
feature  violent  oscillations  of  a  purely  numerical  origin  that  will  eventually  play  havoc  with  the  convergence 
characteristics  of  the  nonlinear  equation  solver. 

Among  the  various  geometric  characteristics  of  shell  equations,  energy  preservation  appears  to  be  the 
most  important  for  the  development  of  robust  time  integration  schemes.  In  fact,  strict  energy  preservation 
at  the  discrete  level  leads  to  unconditional  stability  in  the  nonlinear  regime^  whereas  the  classical  approach 
based  on  the  analysis  of  the  spectral  radius  leads  to  unconditional  stability  in  the  linear  regime  only.  An 
energy  preserving  (EP)  scheme  for  geometrically  exact  shell  is  developed  in  this  paper.  In  addition,  the 
scheme  also  preserves  both  linear  and  anguletr  momenta  of  the  system  at  the  discrete  level.  Unfortunately, 
preservation  of  energy  and  high  frequency  dissipation  cannot  coexist,  imless  energy  is  transferred  from  high 
to  low  frequency  modes,  a  transfer  that  has  no  physical  basis.  Tb  solve  this  problem,  a  family  of  energy 
decaying  (ED)  schemes  that  imply  a  controllable  energy  decay  within  each  time  step  is  proposed  in  this  work. 
In  geometric  ternas,  this  means  that  the  evolution  of  the  system  is  not  confined  to  the  level  set  of  constant 
energy,  but  is  allowed  to  drift  away  from  it  in  a  monotonic  and  controllable  manner.  Since  the  energy  remains 
bounded  at  all  times,  the  scheme  is  unconditionally  stable  for  nonlinear  systems.  Furthermore,  it  can  be 
shown  that  the  energy  dissipation  rtiechanism  of  the  algorithm  is  the  result  of  the  removal  of  the  higher 
frequencies  from  the  computed  response. 

In  related  papers,  various  energy  preserving  and  decaying  geometric  integrators  were  developed  for  rigid 
bodies  and  geometrically  exact  beams  [11,  12,  15,  7,  8],  and  nonlinear  elastodynamics  [10].  The  concepts 
were  extended  to  multibody  systems  featuring  nonlinear  holonomic  constraints  [4,  16,  13,  17,  18,  14];  non- 
holonomic  and  unilateral  constraints  were  treated  in  [6,  5].  The  integration  of  the  present  shell  model  in  a 
general  finite  element  based  multibody  framework  is  discussed  in  [9].  The  proposed  scheme  is  independent 
of  the  choice  of  spatial  discretization  applied  to  the  governing  partial  differential  equations.  In  the  present 
implementation,  the  finite  element  method  is  used,  and  the  mixed  interpolation  of  tensorial  components  [2, 
3,  19]  is  implemented  to  avoid  the  shear  locking  problem.  The  orientation  of  unit  shell  directors  is  described 
by  a  special  family  of  two-parameter  rotations. 

Nonlinearly  energy  decaying  schemes  have  also  been  considered  by  other  authors.  Notably,  a  modification 
of  the  Energy-Momentum  schemes  of  Simo  and  co-authors  was  pursued  for  beams  and  shells  in  refs.  [22,  26], 
following  on  an  idea  proposed  in  ref.  [1]  for  contact  problems.  Unfortunately,  although  not  always  explicitely 
noted  in  the  cited  papers,  the  resulting  algorithms  are  only  first-order  accurate,  and  therefore  of  little  use  in 
practical  applications.  In  refe.  [27,  26],  standard  high  frequency  damping  numerical  integrators  are  combined 
with  the  explicit  enforcement  of  conservation  laws,  in  particular  the  total  energy  with  the  possible  addition 
of  the  system  momenta,  resulting  in  the  Constrained  Energy  Method  (CEM)  or  the  Constrained  Energy- 
Momentum  Algorithm  (CEMA),  respectively.  However,  this  does  not  seem  to  be  an  appropriate  manner 
of  achieving  nonlinear  imconditional  stability:  since  the  higher  modes  are  dissipated  while  the  total  system 
energy  must  remain  constant,  the  algorithm  effectively  transfers  energy  from  the  artificial  high  modes  to  the 
physically  meaningful  low  frequency  modes.  This  un-natural  transfer  was  pointed  out  in  [26]. 

In  contrast  to  the  published  approaches,  the  method  here  proposed  is  the  only  one  to  our  knowledge 
that:  a)  is  higher  order  accurate  (between  two  and  four  [8]),  b)  removes  the  higher  modes,  as  demonstrated 
by  a  spectral  analysis  [8],  while  at  the  same  ensuring  a  strict  algorithmic  total  energy  decay  in  the  nonlinear 
regime. 

The  paper  is  laid  out  as  follows.  The  classical  equations  of  motion  for  geometrically  exact  shells  based 
on  inextensible  unit  directors  are  presented  in  section  2.  Next,  an  EP  scheme  is  developed  in  section  3. 
Section  4  then  pr^ents  an  ED  algorithm  with  tunable  high  frequency  dissipation  that  is  constructed  from 
the  EP  scheme.  Finally,  numerical  examples  are  presented  in  section  5  to  demonstrate  the  efficiency  and 
robustness  of  the  proposed  scheme.  A  discussion  section  concludes  the  paper. 


2 


2  Formulation  of  the  Equations  of  Motion 


2.1  Kinematics  of  the  Shell  Problem 

Consider  a  shell  of  thickness  /i  and  reference  surface  area  fl,  as  depicted  in  fig.  1.  An  inertial  frame  of 
reference  S  consisting  of  three  mutually  orthogonal  unit  vectors  ii,  i2,  is  is  used.  Let  be  the  position 
vector  of  an  arbitrary  point  on  the  reference  surface  of  the  shell,  and  let  C  be  the  material  coordinate  along 
n,  the  normal  to  the  reference  surface.  The  position  vector  r  of  an  arbitrary  point  on  the  shell  in  its  reference 
configuration  is  then 

E(e^^^C)  =  ro(€^a  +  Cn(e^^2)^  (1) 

where  ^^and  are  the  material  coordinates  used  to  represent  the  shell  reference  surface.  The  coordinates 
and  C  form  a  set  of  curvilinear  coordinates  that  are  a  natural  choice  to  represent  the  shell  geoiriCtry. 
The  coordinates  $^and  are  assumed  to  be  lines  of  curvatures  of  the  shell  reference  surface.  The  base 
vectors  are  then 


£2’  Ss 


(2) 


where  Ri  arid  R2  are  the  principal  radu  of  curvature,  =  2:0, o>  the  notation  (•),£,  is  used  to  denote  a 

derivative  with  respect  to  It  is  convenient  to  introduce  a  set  of  three  mutually  orthogonal  unit  vectors 
at  the  shell  reference  surface  (i.e.  at  C  =  0) 

where  Oaa  =  0^ 

Two  fundamental  assumptions  will  be  made  concerning  the  deformation  of  the  shell,  ie.  the  material 
line  initially  normal  to  the  reference  surface  of  the  shell  remains  a  straight  line  and  suffers  no  extension. 
This  is  the  classical  inextensible  director  model.  With  these  assumptions,  th^  position  vector  of  a  material 
point  of  the  shell  writes 

Rie,e,<)=rf,{i\e)+iL{e,e)+cE^{e,e),  (4) 

where  is  the  reference  surface  displacement  vector.  In  the  deformed  configuration,  the  base  vectors 

at  the  shell  reference  surface  are 


da 


G  =  [Gj,  G2,  G3]  =  1,  2,25 

Introducing  the  position  vector,  eq.  (4),  then  yields 


where 


—  [£l)  ^2i  &]  =  £1  + 


..4-^-2  p 

,  £2  -I - ^=^5  ^3 


H- 


^,1  ^,2 

.VSir’  y/^ 


,0 


(5) 


(6) 


(7) 


v^oIT  VS22 

Note  that  is  a  unit  vector,  whereas  and  E.2  are  not  unit  vectors,  nor  are  they  orthogonal  to 

E^,  as  axieil  and  transverse  shearing  strains  develop  during  deformation. 


2.2  Equations  of  Motion 

The  Green-Lagrange  strain  tensor  e  is  defined  as 

e  =  |(G^G-s^p).  (8) 

The  strain  tensor  e  is  defined  in  the  curvilinear  coordinate  S5^tem  defined  by  coordinates  and  G 
However,  it  is  more  convenient  to  work  with  the  strain  tensor  e  defined  in  the  locally  rectanguleir  system 
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defined  by  triad  Cj,  £2,  63,  see  eqs.  (3).  For  shallow  shells  {i.e.  C/^i  <  1  and  C/R2  <  1)  undergoing  large 
displacements  and  rotations  but  small  strains  (all  strain  components  are  assumed  to  be  small  compared  to 
unity),  the  strain-displacement  relationships  can  be  written  as 


e=^[E'^E-I  +  C{E'^H  +  H^E  +  K)], 


(9) 


where 


0  0  ■ 
I/R2  0 
0  0 


(10) 


It  is  clear  that  the  strains  cim  be  expressed  in  terms  of  five  parameters:  the  three  components  of  the 
displacement  field  a  (through  £1  and  £2)  the  two  parameters  defining  the  orientation  of  the  unit 
director  £3.  Virtual  changes  in  the  strain  energy  of  the  structure  are  given  by 


SV  =  f  f  6V  60^1=  f  /  tfe-TdCdfi,  (11) 

Jn  Jh  Jci  Jh 

where  dV  is  the  virtual  strain  energy  density,  and  t  the  second  Piola-Kirchhoff  stress  tensor.  Introducing 
the  strains,  eq.  (9),  and  taking  into  account  the  symmetry  of  the  stress  tensor  then  yields 


6V  =  5E-{E-^(iH)r  +  8H-^ET.  (12) 

The  existence  of  a  strain  energy  density  function  V  is  postulated  here,  hence  the  constitutive  laws  are  of  the 
formT  =  5V/de. 

The  velocity  vector  of  material  point  P  of  the  shell  is  obtained  by  differentiating  the  position  vector, 
eq.  (4),  with  respect  to  time,  to  find  v  =  u  +  The  kinetic  energy  of  the  system  is  now 


/£:  =  jT  jJ^^dCdfi  =  jJ^pjj  udCdn,  (13) 

where  K  is  the  kinetic  energy  density.  Introducing  the  velocity  vector  then  yields 

K  =^p{u-\- C^)-{u  +  C^).  (14) 

Hamilton’s  Principle  can  now  be  expressed  as 


f  ^  /  f  {SK  -  SV)  dCdndt 

Jti  Ja  Jh 

^  ~  lu  lalh  ^  ■{u  +  CE^)  +  5E-{E  +  (,H)t  4-  5H  •  ^Et]  dCdfidt 


=  0.  (15) 


Integrating  through  the  thickness  of  the  shell,  we  get 

{d«  •  [A  -  (iVi,i  +  iV2,2)]  +  5E  ■  [5^  -  (Mi,i  +  M2,2)  +  iV3]  }  dCldt  =  0.  (16) 

In  this  expression,  h  =  mu  +  and  g  =  s*u  -f  are  the  linear  and  angular  momentum  vectors  of 
the  shell,  respectively;  the  mass  coefficients  are  defined  as  m  =  d^,  s*  =  pC  <iC»  I*  =  Sh  <iC-  ‘ 

The  in-plane  forces  are  ^  =  (££^  +  the  out-of-plane  forces  £3  =  EN^,  and  the  bending 

moments  =  {EM^)/^/a^.  The  converted  forces  axe  N*  =  d^,  and  the  convected 

bending  moments  M*  =  [MJ, M2. M3]  = dC- 

The  equations  of  motion  of  shells  could  be  derived  from  this  principle  by  expressing  the  variations  SE^ 
in  terms  of  two  components  of  virtual  rotation. 
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3  Energy  Preserving  Scheme 

Discrete  equations  of  inotion  that  imply  discrete  conservation  laws  for  the  total  mechanical  energy,  linear 
momentum  and  angular  momentum  of  the  system  will  now  be  developed.  The  resulting  scheme  is  closely 
related  to  that  of  ref.  [30],  with  the  slight  difference  of  more  general  definitions  of  the  local  mass  and 
stiffness  tensors  which  are  useful  for  many  engineering  applications.  Similar  energy  preserving  schemes  for 
geometrically  exswjt  shells  with  extensible  directors  are  given  in  ref.  [28]. 

Times  ti  and  tf  denote  the  initial  and  final  times  for  a  time  step,  respectively,  and  the  subscripts  (•)i  and 
(•)/  indicate  quantities  at  tj  and  tf,  respectively.  The  time  step  size  is  denoted  At  =  tf  —  ti.  Furthermore, 
the  subscript  {•)m  is  used  to  denote  mid-point  average  quantities  defined  as 

(■)«.  =  I  !(■)/  + (Oil-  (17) 

The  following  matrix  identity  will  be  used  extensively 

AjBf  -AfBi^  {Af  -  AifBm  +  A^Bf  -  Bi).  (18) 

Hamilton’s  Principle,  eq.  (15),  is  now  approximated  in  time  in  the  the  following  manner 

P  [(M/  -  Mi)  +  C(£3/  “  £31)]  *  ^  ~^^~Af^ 


^i)  ■  (-^m  +  CHm)Ta  +  {Hj  —  Hi)  •  ^EmTa  /  d^dfi  =  0. 


The  change  in  strain  components  from  ti  to  </  is  evaluated  with  the  help  of  identity  (18)  to  find 

e/  -  ei  =  1  [{Ef  -  Ei)^(E„  -b  CHm)  +  {Em  +  CHm)^iEf  -  Ei) 

+{Hf-Hi)^CEm  +  CEl{Hf-Hi)].  (20) 

Over  one  time  step,  the  strain  components  Can  be  approximated  as  e{T})  =  em  +  ri{ef  —  ei)/2,  where  rj  = 
2(t  —  tm)/At  is  the  non-dimensional  time.  If  the  strain  energy  density  function  V  is  viewed  as  a  function  of 
the  scalar  variable  rj,  the  mean  value  theorem  then  implies  the  existence  of  a  ?)  e  [-1,1]  such  that 

This  relationship  defines  the  average  second  Piola-Kirchhoff  stress  tensor,  Tq  =  dV/de\fj.  Combining  this 
result  with  eq.  (20)  then  leads  to 

{Ef  —  Ei)  •  {Em  +  CHm)Ta  +  {Hf  —  Hi)  •  CEmTa  =  (c/  —  Cj)  •  Tq  =  Vf  —  Vi,  (22) 

where  the  symmetry  of  the  stress  tensor  was  taken  into  account.  For  linear  constitutive  laws  of  the  form 
T  =  C*  e,  where  C*  is  the  stiffness  matrix,  the  average  stress  tensor  simply  becomes  Tq  =  C*  Cm- 
The  following  configuration  updates  are  now  defined 


iLf—lLi_  E^f  —  E^i 

At  At 


(23) 


Introducing  eqs.  (22)  and  (23)  into  the  approximate  expression  for  Hamilton’s  Principle,  eq.  (19),  then  leads 
to 


J^f^  {p{ikm  +  CE^rn)-[{Uf-Ui)  +  C{E^f-Esi)]+{Vf-Vi)]  dCdfi 

=  ^  {  f  (*/  +  CEsf)  •  (M/  +  c^f)  - 1  (Mi  +  c4i)  •  (Mi  +  CEsi)  +  {Vf  -  Vi)}  dCdfl 

-^i)  +  iVf-  Vi)]  dCdfi  =  0.  (24) 
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This  result  clearly  implies  the  conservation  of  the  total  mechanical  energy  of  the  system  within  a  step. 
In  sununary,  the  approximate  form  of  Hamilton’s  Principle  given  by  eq.  (19)  leads  to  a  discrete  energy 
conservation  statement,  eq.  (24),  when  the  configuration  updates  are  chosen  according  to  eqs.  (23),  and  the 
average  stress  according  to  eq.  (21). 

Integrating  through  the  thickness  of  the  shell  leads  to 

+  (Msf  ~  Esi)  *  ^ +  M.2m,2)  +  i^mj  |  dfl  =  0.  (25) 

In  this  expression,  the  in-plane  forcra  are  =  {E^NXa  +  the  out-of-plane  forces  = 

and  the  bending  moments  Mam  —  (EmMUa)/ ■  The  discrete  governing  equations  of  motion 
for  shells  are  then 

(26) 

<^l^^-Qm{Mim,l+M2m,2-N^rrd=C^  (27) 

where  p  axe  the  externally  applied  loads,  and  q*  the  externally  applied  moments  measured  in  the  local 
system.  The  finite  change  in  director  oriefitation  E^f  —  E^i  was  expressed  in  terms  of  the  two  parameter 
incremental  rotation  vector,  see  B. 

Invariance  of  the  system  Hamiltonian  under  spatial  translations  and  rotations  implies  the  conservation 
of  the  linear  and  angular  momenta.  Although  discrete  preservation  of  momenta  is  less  crucial  than  discrete 
preservation  of  energy,  it  is  interesting  to  note  that  eqs.  (26)  and  (27)  also  imply  the  discrete  preservation 
of  this  invariant.  At  first,  eqs.  (26)  are  projected  onto  the  test  functions  £  (r^  -l-Um)  and  eqs.  (27)  onto  the 
test  functions  KEsm^  where  tt  is  an  arbitrary  vector  and  i  is  its  associated  skew-symmetric  matrix.  Next, 
integration  over  the  shell  reference  surface  yields 

|i(ro +«m)  •  +i)^2m,2)] 

+  HE^m  *  I  "  ^ (Mlm.l +‘^2m,2) +i^mj  I  dn  =  0.  (28) 
Straightforward  algebraic  manipulations  then  lead  to 

^  ^  [(Eo  +  “/)&/-  (ro  +  li)  hi  +  hmi^f  -  £i) 
where  the  following  result  was  used 

Rirn^Xm  R2mR-2m  EzfimM.2m  d”  EsmE^m  (^®) 

Inserting  the  configuration  updates,  eqs.  (23),  into  eq.  (29)  then  yields 

^  ^  [(To  +  If)  hf  -  (ro  +  Ui)  h,  +  - 13,£.  +  -I-  g^Ezm)]  =  0.  (31) 

It  is  easily  verified  that  hm&m  +£„,^m  —  0-  Hence,  since  x  is  arbitrary,  eq.  (31)  implies  the  discrete 
conservation  of  the  total  angular  momentum,  /^[(^  +  u)  h-f  ££]  dfi.  Finally,  projecting  eqs.  (26)  onto  the 
test  functions  n  and  eqs.  (27)  onto  the  null  test  Actions  gives  the  discrete  conservation  of  the  total  linear 
momentum  /fjh  dfl. 

It  is  important  to  note  that  any  spatial  discretization  of  the  discrete  equations  of  motion  will  inherit  the 
discrete  energy  and  momentum  conservation  statements  just  proved,  when  the  configuration  updates  are 
chosen  according  to  eqs.  (23),  and  the  average  stress  according  to  eq.  (21). 
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4  Energy  Decaying  Scheme 

As  discussed  in  the  introduction,  energy  preservation,  per  5e,  is  not  sufficient  to  yield  robust  time  integration 
schemes.  High  frequency  numerical  dissipation  must  be  added  as  an  inherent  feature  of  the  scheme.  Such  a 
scheme  will  now  be  constructed  for  the  shell  equations  of  motion  using  the  EP  scheme  as  a  basic  building 
block.  Further  details  on  the  scheme  are  given  in  ref.  [8]. 

First,  an  additional  state  is  introduced  at  time  tj  =  lime_»o(^i  +  €),  and  the  subscript  (•)j  is  used  to 
denote  quantities  at  this  time.  The  following  averages  are  now  defined 

(•)s  =  |lO/  +  Oi];  (•)/.  =  I  ((•),  + (Oil .  (32) 

The  ED  scheme  proceeds  from  the  initial  to  the  final  time  by  means  of  two  coupled  steps:  one  step  from  U 
to  t/,  the  other  from  U  to  tj.  The  time-discrete  equations  of  dynamic  equilibrium  are 


hf-hi 


-  {E.\g,l-^R2g,2)  =Pg> 


qt  zi _ ^  _  qT 


+  ;  [(&,.!  +liw)  -  (w.,,. +&,.2)1  -  a; 


At  ^  3 


+  i  [Ql  (M,,,,  +  -  N^g)  -  Ql  (Mip.i  +  M2p,2  -  iVgp)]  =  s:- 


The  configuration  update  relationships  are  given  as 

My  ==  2ii  +  At  -f  —  [uf  —  -  a{Uj  —  Ui)] /Q; 

r  .  .  .  *1  (^5) 

£3/  =  +  At  (^y  -f  Msj)/2,  E^j  =  -  At  -  Oi{E^j  ~  E3i)\ 

where  a  is  a  tuning  parameter  that  controls  the  amount  of  numerical  dissipation  provided  by  the  scheme, 
while  the  forces  and  moments  are  given  by 

Nap  =  Kah+<^  {Kaj  ^  Nai)/^\  Map  =  Ma/.  +  ^  (Ma,*  ^  Mai)/2.  (36) 

Using  developments  similar  to  those  exposed  for  the  EP  scheme,  it  can  be  easily  shown  that  the  proposed 
discrete  equations  imply 

iKf  +  Vf)-{Ki  +  Vi)  +  a(^  =  0.  (37) 

c®  is  a  positive  quantity  given  by 

=  l^l[rn\\u\\-\\u\\+2s*\\u\\-\\^\\+rm\\-m\\]dQ 
+  U  ||e||C*|U||dn>0,  (38) 


where  ||  •  ||=  (•)y  —  (•)i  is  the  jump  between  U  and  tj.  This  result  implies  the  decay  of  the  total  mechanical 
energy  over  one  step  of  the  algorithm,  {Kf  +  Vf)  <  {Ki  -I-  Vi).  The  parameter  a  clearly  controls  the  amount 
of  energy  that  is  dissipated  within  the  step.  Two  such  parameters  could  be  used,  controlling  the  amount  of 
dissipated  kinetic  and  strain  energies,  respectively,  but  this  level  of  complexity  does  not  seem  to  be  necessary. 
The  property  of  preservation  of  momentum  observed  in  the  EP  case  is  lost  in  the  ED  algorithm. 

If  the  above  ED  scheme  is  applied  to  a  single  degree  of  freedom  linear  oscillator,  the  asymptotic  value 
of  the  spectral  radius  of  the  amplification  matrix,  pc»,  is  found  to  be  poo  =  (1  ck)/(1  +  a).  For  a  =  1, 

Poo  =  0,  and  asymptotic  annihilation  is  achieved.  If  a  =  0,  Poo  =  I,  and  in  view  of  eq.  (37),  energ>^  is  exactly 
preserved.  Hence,  the  ED  scheme  is  in  fact  a  family  of  schemes  with  a  single  tuning  parameter,  a,  that 
controls  the  amount  of  high  frequency  numerical  dissipation;  both  asymptotic  annihilation  or  exact  energy 
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preservation  can  be  achieved  with  the  same  scheme  by  using  a  =  1  or  0,  respectively.  The  spectral  radius, 
algorithmic  damping  ratio  and  relative  period  error,  as  defined  in  ref.  [24],  are  plotted  in  figs.  2  through  4 
for  varying  q. 

An  error  analysis  based  on  Taylor  series  expansions  for  the  linear  harmonic  oscillator  model  problem 
determined  that  the  ED(a  =  0)  scheme  is  fourth  order  accurate.  For  the  same  model  problem,  ED(a  =  1)  is 
third  order  accurate.  For  general  nonlinear  problems,  numerical  experiments  show  that  all  methods  exhibit 
convergence  rates  ranging  between  two  and  three. 

5  Numerical  Examples 

All  the  examples  described  in  this  section  will  be  treated  with  the  proposed  ED  family  of  schemes  corre¬ 
sponding  to  values  of  the  tuning  parameter  a  €  [0, 1).  Although  any  value  of  a  within  this  range  can  be 
used,  the  examples  described  here  will  contrast  the  two  extreme  choices.  For  a  =  1  (poo  =  0),  asymptotic 
annihilation  is  obtained.  On  the  other  hand,  for  a  =  0  (poo  =  1)?  exact  energy  preservation  is  achieved. 

5.1  Clamped  Half-Cylinder  under  Point  Load 

Consider  the  half-cylinder  of  radius  ii  =  1.2  m  and  width  6  =  2  m  depicted  in  fig.  5.  The  shell  has  a 
thickness  t  =  6  mm,  is  build-in  along  edge  BC  and  free  along  the  other.  The  structure  is  made  of  aluminum; 
Young’s  modulus  E  =  73  GPa,  Poisson’s  ratio  i/  =  0.30  and  density  p  ==  2700  kg/m^.  At  point  D,  the  shell 
is  subjected  to  a  concentrated  load  P  —  ~Po(t)(ii  + 12+13)-  The  magnitude  of  the  load  is 

P0(t)  =  | 

where  P  =  0.1  kN  and  T  =  2.0  s.  The  shell  was  modeled  by  a  regular  8x4  mesh  of  quadratic  elements.  All 
simulations  were  run  with  a  time  step  At  =  5.0  10~®^  s,  for  a  total  time  of  6  s. 

Under  the  effect  of  the  applied  loads,  the  shell  bends  predominantly  in  the  vertical  direction,  its  direction 
of  least  bending  stiffness,  as  illustrated  in  fig.  6  that  shows  the  configuration  of  the  system  at  various  instants 
in  time.  The  three  components  of  displacements  at  point  D  are  shown  in  fig.  7;  vertical  displacements  of  up 
to  0.6  m  are  observed.  The  components  of  transverse  shearing  forces  axe  shown  in  fig.  8. 

Next,  the  same  problem  was  simulated  with  poo  —  1»  i-c.  with  no  high  frequency  dissipation.  The 
corresponding  results  are  shown  in  figs.  7  to  8.  Displacement  and  moment  results  are  found  to  be  in 
excellent  agreement.  At  the  scale  of  the  figures,  they  are,  in  fact,  indistinguishable.  For  the  period  of  time 
2  <  t  <  6  s,  the  system  is  not  subjected  to  any  loading,  and  the  total  mechanical  energy  of  the  system 
should  remain  constant.  For  the  ED  (a  =  0)  scheme,  the  energy  is  indeed  preserved,  as  expected;  for  the 
ED(a  =  1)  scheme,  0.3%  of  the  energy  is  numerically  dissipated  in  this  period,  as  depicted  in  fig.  9.  It 
could  be  concluded  that  the  ED(a  =  0)  and  ED(q  —  1)  solutions  are  nearly  identical,  and  that  numerical 
dissipation  is  not  necessary.  However,  the  ED(a  =  1)  and  ED(a  =  0)  scheme  predictions  for  the  transverse 
shearing  forces,  shown  in  fig.  8,  are  markedly  different.  ED(a  =  0)  predictions  for  force  components  P13 
show  high  frequency  oscillations  that  are  absent  in  the  corresponding  ED(a  =  1)  predictions.  A  simulation 
using  the  ED(a  =  1)  scheme  with  a  time  step  At  ^  1.0  10“^®^  s  showed  that  the  ED(o[:  =  1)  predictions 
are  converged.  A  simulation  using  the  ED(a  =  0)  scheme  and  the  same  smaller  time  step  yielded  results 
with  increased  high  frequency  oscillations  for  the  force  predictions.  It  should  be  noted  that  the  dynamic 
response  of  this  simple  system  is  very  smooth;  yet  even  here,  high  frequency  numerical  dissipation  appears 
to  be  necessary  to  obtain  a  smooth,  converged  solution. 

6.2  Dynamic  Response  of  a  Plate  with  Edge  Beams 

Consider  the  rectangular  plate  of  length  L  2  m,  width  6  =  0.05  m  and  thickness  t  =  2  mm  as  depicted  in 
fig.  10.  Two  circular  beams  of  radius  r  =  2  mm  are  attached  at  the  plate  edges.  A  third  circular  beam  is 
located  at  the  center  of  the  plate.  All  components  are  made  of  aluminum;  Young’s  modulus  E  =  73  GPa, 
density  p  —  2700  kg/m^.  The  total  mass  of  the  edge  beams  is  10  kg  each,  and  that  of  the  central  beam  is 
1  kg.  The  plate  is  subjected  to  uniformly  distributed  loads  Fm  and  Fi  along  FE  and  CJ5,  respectively.  The 
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components  of  these  loads  along  the  ij  and  23  axes  are  Fim  =  40  N/m  and  F^m  —  80  N/m,  respectively,  and 
Fu  =  “20  N/m  and  Fsi  =  —60  N/m,  respectively.  The  common  time  history  of  each  loading  component  is 


(1  “  cos27rt/r)/2  t  <  r, 

t  >r, 


where  T  =  3.0  s.  The  plate  is  discretized  with  10  quadratic  plate  elements  along  its  length.  A  constant  time 
step  At  =  6  10"®^  s  was  used  for  all  simulations. 

At  time  t  >  3  s  the  applied  load  vanishes,  and  the  system  total  mechanical  energy  should  remain  constant. 
The  evolution  of  the  energy  of  the  system  for  three  different  cases  is  shown  in  fig.  11:  the  plate  model  using  the 
ED(a  =  1)  scheme,  the  same  plate  model  using  the  ED(q  =  0)  scheme,  and  a  simplified  model  of  the  system 
using  beam  elements  and  the  ED(a  =  1)  scheme.  For  times  t  >  3  s,  the  total  energy  remains  a  constant  for  the 
ED(q  =  0)  scheme  and  nearly  constant  for  the  ED(q:  ~  1)  schemes.  The  relative  energy  loss  is  also  presented 
in  the  figure:  0.6%  of  the  energy  was  dissipated  by  the  ED(a  =  1)  scheme  in  the  period  t  e  [3, 20]  s.  This 
figure  clearly  demonstrates  the  non-increasing  property  of  the  energy  evolution  for  the  proposed  ED  (a  =  1) 
schemes,  as  opposed  the  constant  energy  predicted  in  ED(a  =  0)  simulations.  The  trajectory  of  the  plate 
mid-span  point  is  shown  in  fig.  12:  good  correlation  is  observed  between  the  predictions  of  the  three  models. 
The  beam  model  is  slightly  off  due  to  the  inherent  simplifying  assumptions.  The  behavior  of  the  quarter-span 
axial  force  and  transverse  shear  force  are  shown  in  fig.  13  and  14,  respectively.  The  poor  predictions  of  the 
ED(a  =  0)  schemes  are  obvious  in  these  two  plots.  The  history  of  axial  force  presents  violent  oscillations 
with  amplitudes  an  order  of  magnitude  larger  than  those  observed  for  the  ED  (a  =  1)  scheme.  The  history 
of  the  transverse  shear  force  predicted  by  the  ED(a  =  0)  scheme  quickly  diverges  from  the  ED(a  =  1) 
predictions  for  both  beam  and  plate  models.  To  ascertain  the  accuracy  of  the  ED(a  =  1)  predictions,  a 
convergence  study  was  performed.  Nearly  identical  results  were  found  with  smaller  time  step  sizes  At  =  3.0 
and  1.0  lO"”®^  s,  or  when  using  the  time  adaptivity  procedure.  On  the  other  hand,  oscillations  of  increasing 
amplitude  were  found  as  the  time  step  size  is  reduced  in  the  ED  (a  =  0)  scheme.  Furthermore,  the  time 
adaptivity  procedure  failed  to  yield  any  results  because  the  time  step  size  was  driven  to  unreasonably  small 
values,  At  =  s,  as  the  procedure  tries  to  cope  with  increasingly  violent  oscillations. 


5.3  Dynamic  Response  of  a  Cruciform 

Consider  a  cruciform  consisting  of  four  thin  panels  {Panels  A,  5,  C,  and  D)  connected  to  a  central  beam,  as 
depicted  in  fig.  15.  Each  panel  is  of  thickness  t  =  4  mm,  length  L  =  1.2  m,  and  width  6  =  0.1  m.  The  central 
beam  has  a  square  cross-section  of  width  a  =  8  mm.  A  mass  M  =  12  kg  is  attached  at  the  tip  of  the  central 
beam  at  point  T.  Panels  and  beam  are  simply  supported  at  the  root  of  the  cruciform.  A  concentrated  load 
P(t)  is  applied  at  point  T.  The  load  acts  in  the  plane  defined  by  axes  12  and  and  makes  a  30  degree  angle 
with  axis  12-  All  components  are  made  of  aluminum  with  properties  given  in  the  previous  example.  The 
time  history  of  the  applied  load  is 


p(<)  =  | 


Po  (1  -  cos  2irtlT)l2  t<T, 
0  t>T, 


where  Po  =  1.2  kN  and  T  =  0.1  s. 

As  the  applied  load  increases,  in-plane  stresses  in  the  panels  rapidly  increase  and  buckling  takes  place 
in  those  panels  subjected  to  compression,  as  can  be  observed  in  fig.  16  that  depicts  the  configuration  of 
the  cruciform  at  two  instants  in  time.  The  trajectory  of  point  T  projected  onto  plane  12,  is  is  shown  in 
fig.  17.  For  reference,  the  corresponding  trajectory  of  a  beam  with  cross-sectional  properties  equivalent  to 
those  of  the  cruciform  is  also  presented.  Of  course,  the  equivalent  beam  model  is  much  stiffer  since  it  does 
not  allow  buckling  to  take  place.  Furthermore,  the  motion  remains  confined  to  the  plane  defined  by  axis  ij 
and  the  line  of  action  of  the  applied  load.  When  each  panel  is  modeled  individually,  the  stiffness  of  system 
varies  both  spatially  and  temporally,  giving  rise  to  the  more  complex  motion  shown  in  fig.  17.  The  total 
mechanical  energy  of  the  system  is  shown  in  fig.  18.  Prom  time  t  =  0.1  to  0.2  s,  the  system  is  free  and  its 
total  mechanical  energy  should  remain  constant.  Due  the  dissipative  nature  Of  the  integration  scheme,  a 
small  amount  of  energy  is  dissipated  over  that  period  of  time:  2.7%  of  the  energy  was  dissipated  over  the 
2435  time  step  period. 


The  root  shear  force  and  quarter-span  bending  moment  in  Panels  A  and  C  are  shown  in  fig.  19  and  20, 
respectively.  Each  panel  undergoes  alternating  phases  of  tensile  and  compressive  loading.  During  the  com¬ 
pressive  phases,  buckling  takes  place,  and  large  shear  forces  and  bending  moments  are  observed  in  contrast 
with  the  tensile  phases  during  which  these  quantities  remain  much  smaller. 

5.4  Snap-Through  of  a  Cylindrical  Shell 

The  snap-through  behavior  of  a  cylindrical  shell  under  a  concentrated  load  was  investigated  in  ref  [27].  The 
shell  consists  of  a  60  degree  sector  of  a  cylinder  of  height  h  =  5  m,  radius  i?  =  5  m  and  thickness  t  =  0.1  m, 
as  shown  in  fig.  21.  Material  properties  are:  Young’s  modulus  E  —  210  GPa,  Poisson  ratio  u  =  0.25  and 
density  p  =  10^  kg/m^.  The  two  straight  edges  of  the  shell  are  simply  supported,  while  the  two  curved  edges 
are  free. 

A  concentrated  force  F  is  applied  at  the  shell’s  apex.  This  force  linearly  increases  from  0  to  5  10*^  N  in 
0.2  s,  then  is  held  constant  at  that  value.  The  simulation  ends  at  time  t  =  0.3  s.  Due  to  the  symmetry  of 
the  problem,  a  quarter  shell  only  is  modeled;  a  regular  4x4  mesh  of  quadratic  elements  was  used.  The  time 
step  size  was  selected  as  At  =  10"^  s. 

As  the  load  increases,  the  shell  apex  displacement  increases,  then  suddenly,  snap-through  takes  place  and 
curvature  reverses.  Curvature  reversal  initiates  in  the  region  of  the  applied  load,  then  quickly  propagates 
throughout  the  entire  structure,  which  undergoes  subsequent  violent  oscillations.  Snapshots  of  the  system 
at  various  instants  in  time  are  given  in  fig.  22.  The  vertical  displacements  of  the  point  of  application  of 
the  load  computed  with  the  ED(a  =  1)  scheme  is  shown  in  fig.  23.  Note  the  gradual  increase  of  the  shell 
deflection,  until  collapse  at  buckling  and  the  resulting  vibratory  response  in  the  inverted  configuration. 

Ref.  [27]  presents  simulations  of  this  problem  using  various  schemes:  the  generalized-Q  [21]  and  the 
CEMA  [27]  schemes.  The  former  scheme  features  high  frequency  numerical  dissipation  and  linear  stability 
properties,  while  the  latter  adds  to  the  genera;lized-a  method  a  constraint  on  the  total  mechanical  energy  of 
the  system.  CEMA  is  therefore  both  energy  preserving  and  high  frequency  dissipative.  The  results  presented 
in  fig.  23  are  in  close  agreement  with  those  obtained  with  the  generalized-a  scheme,  but  quite  different  from 
those  predicted  by  CEMA  in  the  post  buckling  regime.  It  is  important  to  realize  that  the  higher  modes  are 
only  an  artifact  of  the  discretization  process,  and  should  therefore  be  removed  from  the  computed  response. 
A  standard  scheme  like  the  generalized-a  method  accomplishes  this  goal  through  the  characteristic  low-pass 
shape  of  its  spectral  radius;  however,  there  is  no  guarantee  that  energy  will  not  be  allowed  to  grow  within  one 
step  for  nonlinear  problems.  In  contrast,  CEMA  enforces  the  exact  conservation  of  energy  in  the  nonlinear 
regime,  but  at  the  same  time  inherits  high  frequency  dissipation  from  the  underlying  generalized-a  algorithm. 
Consequently,  an  artificial  mechanism  for  transferring  energy  from  the  higher  (artificial)  modes  to  the  lower 
modes  is  created  that  drives  the  response  to  an  erroneous  solution.  In  contrast,  the  proposed  ED  scheme 
achieves  both  nonlinear  stability  and  high  frequency  dissipation. 

Next,  the  same  problem  was  simulated  with  poo  =  with  no  high  frequency  dissipation.  In  this 

case,  two  refinements  in  time  step  size  were  required  to  successfully  complete  the  simulation,  one  at  time 
t  =  0.1665  s  (At  =  5  10”^  s),  the  other  at  time  t  =  0.2142  s  (At  =  2.5  10“^  s).  Deflections  predicted  by  the 
ED(a  =  0)  and  ED(a  =  1)  schemes,  shown  in  fig.  23,  are  in  good  agreement  during  the  initial  snap-through 
phase,  but  become  increasingly  different  during  the  subsequent  oscillations.  The  force  and  velocity  fields  are 
markedly  different.  The  plate  center  forces  for  both  ED(a  =  0)  and  ED(a  =  1)  schemes  are  shown  in  fig.  24. 
The  forces  predicted  by  the  ED(a  =  0)  scheme  present  violent  oscillations  of  amplitude  up  to  an  order  of 
magnitude  larger  than  those  predicted  by  the  ED(a  =  1)  scheme.  These  violent  oscillations  hamper  the 
convergence  of  the  Newton  process  at  each  time  step,  leading  to  the  need  for  smaller  time  steps.  The  same 
observations  can  be  made  about  fig.  25  which  compares  the  plate  center  vertical  velocity.  Violent  oscillations 
are  initiated  at  snap-through  and  the  strict  preservation  of  energy  implied  by  the  ED(a  =  0)  scheme  prevents 
any  subsequent  decay  of  these  vibrations.  Since  vibratory  stresses  are  a  great  importance  to  designers,  it  is 
essential  to  assess  the  ability  of  new  integration  schemes  to  reliably  predict  these  quantities.  It  is  unfortunate 
that  many  scientific  publications  about  geometric  integration  only  present  responses  for  preserved  quantities 
such  as  total  mechanical  energy  or  momentum.  The  above  plots  demonstrate  that  while  ED  (a  =  0)  scheme 
might  perform  very  well  for  the  prediction  of  total  energy,  momentum,  or  even  displacement  fields,  they  are 
unable  to  reliably  predict  other  important  fields  such  as  velocities  and  internal  stresses.  Consequently,  such 
schemes  are  of  little  values  in  real  life  applications. 
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6  Conclusions 

In  this  worlc,  a  robust  algorithm  for  the  dynamic  analysis  of  geometrically  exact  shell  structures  was  presented. 
The  method  is  geometry-based,  i.e.  it  incorporates  knowledge  about  specific  qualitative  features  of  the 
underlying  partial  differential  equations.  However,  departing  from  the  classical  approaches  based  on  strict 
preservation  of  energy,  the  method  presented  here  allows  the  system  to  drift  away  from  the  level  set  of 
constant  energy  in  a  controlled  and  tunable  manner. 

This  feature  achieves  two  goals.  First,  a  bound  is  placed  on  the  total  mechanical  energy  of  the  discrete 
system,  leading  to  the  concept  of  nonlinear  unconditional  stability;  this  stability  criterion  is  stronger  than 
that  obtained  through  the  classical  analysis  of  numerical  schemes.  The  resulting  numerical  procedure  is 
endowed  with  superior  robustness,  an  important  feature  when  dealing  with  complex  engineering  problems. 
Second,  the  monotonic  energy  drift  is  associated  with  numerical  dissipation  of  the  high  frequency  modes. 
This  tunable  dissipation  makes  the  algorithm  stiffly  accurate,  and  avoids  the  build  up  of  energy  in  the  higher 
modes  that  are  an  artifact  of  the  spatial  discretization  process. 

The  proposed  scheme  can  deal  with  general  shell  structures  and  is  not  tied  to  a  specific  spatial  dis¬ 
cretization  of  the  governing  partial  differential  equations.  Kinematic  nonlinearities  are  treated  in  a  rigorous 
manner,  and  material  nonlinearities  can  be  handled  when  the  constitutive  laws  stem  from  the  existence  of  a 
strain  energy  density  function.  The  efficiency  and  robustness  of  the  proposed  approach  were  demonstrated 
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A  Rodrigues  Parameters 

It  is  clear  that  any  parameterization  of  finite  rotations  [25]  could  be  used  with  the  present  shell  formulation. 

finifp  rT!"’  Rodrigues  parameters  r  =  2k  tan<i/2,  where  is  the  magnitude  of  the 

famte  rotation  and  k  the  components  of  the  unit  vector  about  which  it  takes  place.  The  following  notation 
IS  mtroduced  r„  =c^^/2^1/  (1  r/4),  and  the  fimte  station  tenter  fi  then  write 

R{l)  =  /  +  To  r  +  y  rr.  (1^ 

The  following  decomposition  of  the  rotation  tensor  is  extensively  used  in  this  work 

B  Orientation  of  a  Unit  Director 

Consider  a  unit  vector  ig,  cafied  a  director,  that  rotates  to  a  final  orientation  e,.  For  convenience  this 
director  is  considered  to  be  the  third  unit  vector  of  a  triad  S  defined  by  ij,  ig,  rotating  to  a  triad  5*’ with 

rotation  t”  ■*’  two  triads  is  where  i?  is  an  orthogonal 

r  t  director,  this  rotation  tensor  is  not  uniquely  defined,  as  any 

rotation  about  the  director  leaves  its  orientation  unchanged.  A  virtual  change  in  the  director  orientation  is 

5e3  =  ef^, 

where  6ip  is  the  virtual  rotation  vector,  ^  =  SRRT. 

The  components  of  the  virtual  change  in  director  orientation  measured  in  S*  become 


R  Se^  =  R'^'^Sip  —  ’i^R'^S'ijj  =  't^6ij}*  =  Srjii^  ^  /2i 

~  0 

components  of  the  virtual  rotation  vector  in  5*.  This  relationship  clearly  demonstrates 

liH  to  virtual  rotations  of  the  director  about  its  own  orientation, 

will  not  affect  virtual  changes  m  the  director  orientation,  and  hence,  setting  =  0  is  a  valid  choice  The 
following  notation  is  adopted  -  a  -Ps  ^  vcinu  cnoice.  me 

^  ^  =  [ill  12]*  (3) 

Sa  IS  a  2  X  1,  “two  parameter”  virtual  rotation  vector.  It  follows  that  6^  =  R  Sip*  =  Rb  6a*,  and  hence 

Se^=R7^bSa*.  (4) 

If  Rodrip^  parameters  are  used  to  parameterize  R,  an  equivalent  expression  can  be  obtained  for  finite 
changes  m  director  orientation  with  the  help  of  eq.  (2) 

Sfl./ ~  =  Rm  ffl,  6  a  ~QmS*’,  £*  =  6s*, 

where  r*  are  the  Rodrigues  parameter  measured  in  5*,  ands*  the  corresponding  “two  parameter”  incremental 
rouation  vector. 
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Figure  2:  Spectral  radius  of  the  ED  scheme  for  varying  a:  o  =  0.0,  □  =  0.2, 0  =  0.4,  A=  0.6,  V  =  0.8,  >=1.0, 
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Figure  3:  Algorithmic  damping  ratio  of  the  ED  scheme  for  varying  a:  o  =  0.0,  □  =  0.2, 0  =  0.4,  A=  0.6,  V 
0.8,  >=1.0. 
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Figure  6:  Configuration  of  the  system  at  various  instants  in  time. 


Figure  7:  Displacement  components  at  point  D,  Ui:  A;  U2*  o;  U^:  □.  ED(a  =  1)  scheme:  solid  line; 
ED(a  =  0)  scheme:  dashed  line. 


I 
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Figures:  Time  history  of  transverse  shear  forces  at  point  M.  F13:  top  figure;  F23:  bottom  figure.  ED(q;  =  1) 
scheme:  solid  line;  ED(a  =  0)  scheme:  dashed  line. 
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Figure  9:  Time  history  of  total  energy  of  the  system  (top  figure).  Relative  total  energy  loss  from  time 
t  e  [2,6]  sec  (bottom  figure).  ED(a  =  I)  scheme:  solid  line;  ED(a  =  0)  scheme:  dashed  line. 
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TIME  [sec] 


Figure  11:  Time  history  of  the  total  mechanical  energy  (top  figure)  and  relative  energy  loss  (bottom  figure). 
Plate  model  (ED(a  =  1)):  solid  line;  Plate  model  (ED(a  =  0)):  dashed  line;  Beam  model  (ED(a  =  1)): 
dashed-dot  line. 
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Figure  13:  Time  history  of  the  quarter-span  axial  force.  Plate  model  (ED(a  =  1)):  solid  line;  Plate  model 
(ED(a  =  0)):  dashed  line;  Beam  model  (ED(a  =  1)):  dashed-dot  line. 
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Figure  14:  Time  history  of  the  quarter-span  transverse  shear  force.  Plate  model  (ED(a  =  1)):  solid  line: 
Plate  model  (ED(a  =  0)):  dashed  line;  Beam  model  (ED(a  =  1)):  dashed-dot  line. 
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Figure  15:  Configuration  of  the  cruciform  problem. 
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Figure  17:  Trajectory  of  point  T  of  the  cruciform  projected  in  the  plane  jjt  is-  Solid  line:  shell  model; 
dashed  line:  beam  model. 
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Figure  18:  Total  mechanical  energy  of  the  system.  Solid  line:  shell  model;  dashed  line:  beam  model. 
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ROOT  TRANSVERSE 


Figure  19:  Time  history  of  the  root  shear  force  in  Panel  A  (solid  line)  and  Panel  C  (dashed  line). 
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Figure  21:  Configuration  of  the  snap-through  problem. 
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[  t  =  0.06  IsecJ 


^  t  =  0.12  [sec] 
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Figure  23:  Time  history  of  the  plate  center  vertical  displacement.  ED(a  =  1)  scheme:  solid  line;  ED(a  =  0) 
scheme:  dashed  line. 
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Figure  24:  Time  history  of  the  plate  center  forces  for  the  ED(a  =  1)  and  ED(a  =  0)  schemes.  Force  Fi: 
solid  line;  F3:  dashed  line.  For  clarity  the  ED(a  =  0)  results  are  shifted  downwards  3  10®  N. 
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